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Preface 


A 


The  problem  of  modeling  optical  turbulence  effects  within  target  acquisition 
models  has  largely  been  limited  to  characterizing  the  influence  of  turbulent 
blur  on  the  modulation  transfer  function  (MTF)  for  the  system  under  study. 
However,  as  useful  as  this  characterization  may  be,  the  effects  of  turbulence  are 
often  not  limited  to  the  MTF  alone.  The  MTF  is  only  one  of  three  turbulence 
effects  on  impacting  sensing  systems. 

The  remaining  two  effects  are  normally  described  as  scintillation  and  angle-of- 
arrival  variations.  Scintillation  causes  temporal  fluctuations  in  received  intensity 
across  a  target  and/or  background.  The  most  obvious  example  of  scintillation 
can  be  seen  when  approaching  an  oncoming  vehicle  on  the  highway.  A  glint 
portion  of  the  vehicle  (usually  a  sun  reflection  from  the  vehicle  windshield) 
appears  to  fluctuate  in  brightness.  Scintillation  is  usually  unimportant  for 
objects  of  near-uniform  brightness,  but  the  appearance  of  hot  spots  (glint 
features)  will  be  significantly  influenced  by  scintillation  effects.  The  influence 
of  scintillation  can  be  considered  a  contributor  to  the  overall  noise  within  an 
atmosphere-optics-electronics-human  detection  system. 

Further,  system  designers  tend  to  develop  systems  which  suppress  scintillation 
effects.  The  log  variance  of  scintillation  varies  as  A;7/6,  where  k,  is  the  radiation 
wavenumber  27r/A.  Systems  operating  at  longer  wavelengths  experienced 
significantly  reduced  scintillation  effects  compared  to  visible  systems.  Receiver 
aperture  diameter  also  influences  the  scintillation.  For  large  receiver  optics  the 
scintillation  is  aperture  averaged.  Hence,  more  scintillation  is  seen  by  the  human 
eye  than  through  a  telescope  at  the  same  wavelength.  There  is,  however,  a 
tradeoff  between  scintillation  and  atmospheric  blurring  effects.  Larger  apertures 
that  tend  to  reduce  scintillation  also  tend  to  cause  increased  turbulent  blurring. 

At  infrared  (IR)  wavelengths  both  scintillation  and  blurring  effects  are  reduced, 
yet  turbulence  effects  are  still  a  factor  due  to  image  distortion.  In  the  far 
IR,  angle-of-arrival  fluctuations  are  the  more  pervasive  turbulence  influence. 
Also  known  as  image  wander,  shimmer,  heat  boil,  and  jitter,  angle-of-arrival 
fluctuations  are  the  first  noticable  turbulent  distortion. 

Currently,  angle-of-arrival  influences  are  not  considered  in  target  acquisition 
models.  The  reason  for  this  deficiency  is  rather  simple:  Systems  analysis  models 
generally  utilize  linear  shift  invariant  (LSI)  assumptions.  Angle  of  arrival  effects 
are  not  LSI  because  of  a  related  parameter  called  the  isoplanatic  patch  size, 
0O.  When  6q  is  large,  scene  objects  appear  to  move  in  position  due  to  angle-of- 
arrival  fluctuations,  yet  tend  to  maintain  their  structural  coherence  because  the 
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same  error  applies  to  every  point  in  the  object.  But  as  do  shrinks  in  angular 
extent  below  the  angular  subtense  of  an  object  of  interest,  different  portions 
of  the  object  appear  to  shift  independently  in  angle  relative  to  one  another. 
The  effects  of  such  a  distortion  on  an  observer’s  ability  to  percieve  the  object 
then  involves  a  complex  interaction  between  object  characteristics,  background 
characteristics,  and  the  qualities  of  the  perturbations.  It  becomes  necessary 
to  conduct  observer  tests  to  determine  the  gestalt  of  the  entire  process  on  the 
ability  of  the  observer  to  understand  the  fluctuating  scene.  In  the  end  it  may 
be  possible  to  assign  complexity  values,  similar  to  the  o  values  used  in  target 
acquisition  codes,  to  evaluate  the  degrading  effects  of  the  image  distortion  on 
the  observer  capabilities.  These  observer  tests  will  require  simulation  of  image 
distortion  through  an  image  modification  code. 

This  paper  provides  support  for  calculations  needed  in  such  image  distortion 
codes.  In  it,  we  describe  a  method  of  modifying  images  for  angle-of-arrival 
distortion  and  for  the  more  standard  phase  screen  approach.  To  perform 
either  of  these  tasks  properly  requires  knowledge  of  the  appropriate  spatial 
frequency  spectrum  of  the  turbulent  fluctuations,  and  knowledge  of  the  means 
of  representing  these  fluctuations  in  a  model  of  the  atmosphere.  One  must 
then  define  a  method  of  image  manipulation  which  supports  the  calculation  of 
distorted  images  from  undistorted  source  images.  This  text  discusses  use  of 
a  turbulence  spectrum  developed  recently  which  handles  both  inner  and  outer 
scale  influences  on  the  refractive  index  spectrum  within  the  earth’s  surface  layer 
atmosphere. 
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Executive  Summary 


Rationale 

Forward-Looking  Infrared  (FLTR)  devices  have  been  a  part  of  the  Army 
inventory  for  over  twenty  years.  In  that  time  the  state-of-the-art  in  FLIR 
technology  has  improved  significantly.  We  arc  now  using  third-generation  FLIR 
devices,  which  exhibit  significantly  lower  noise  and  higher  spatial  resolution  than 
was  possible  twenty  years  ago.  However,  with  that  increased  resolution  has  come 
the  realization  that  turbulence  effects  are  not  minimal  when  compared  with  the 
system  noise  in  these  new  sensors.  Further,  the  next  generation  of  FLIRs  will 
likely  involve  higher  magnification  devices  designed  for  specialty  applications 
such  as  narrow  fields  of  regard.  To  enable  such  applications,  the  sensor  aperture 
diameter  may  need  to  be  increased,  leading  to  larger  amounts  of  turbulent  blur 
on  the  received  imagery.  Lower  noise  thresholds  and  higher  resolution  detection 
grids  often  means  being  able  to  detect  image  distortion  due  to  angle-of-arrival 
fluctuations. 

Because  of  these  concerns,  and  because  current  sensor  performance  prediction 
models  do  not  incorporate  turbulence  effects,  it  is  reasonable  to  attempt  to 
simulate  turbulence  effects  in  a  modeling  environment.  Traditionally,  analysis 
of  turbulence  effects  has  followed  the  course  of  incorporating  turbulent  blur 
effects  as  a  multiplicative  atmospheric  turbulence  Modulation  Transfer  Function 
(MTF)  when  combined  with  the  system  MTF.  However,  for  FLIR.  systems  the 
chief  turbulence  concern  may  be  apparent  object  shape  distortion  and  this  effect, 
because  it  is  not  shift  invariant,  cannot  be  modeled  via  an  MTF.  Rather,  it  must 
be  modeled  using  image  processing  techniques. 


Analyses 

To  support  these  image  processing  techniques,  one  must  develop  a  methodology 
akin  to  methods  originally  developed  for  high  energy  laser  simulations.  The 
principle  method  developed  in  the  1970’s  to  study  these  propagation  problems 
involved  mathematical  constructs  known  as  phase  screens.  While  the  literature 
on  phase  screens  is  somewhat  extensive,  certain  key  elements  of  the  analyses 
in  this  literature  involve  undeveloped  assumptions,  certain  mathematical  flaws, 
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and  exercises  left  to  the  reader.  This  approach  did  not  seem  satisfactory,  leading 
to  the  development  of  this  document. 

In  this  document  we  explore  the  basic  equations  necessary  to  describe 
Fourier  transforms,  Fourier  series  expansions,  and  the  fundamental  equations 
undergirding  the  Fast  Fourier  Transform  technique.  Because  these  methods 
are  essentially  related,  we  show  the  means  of  transforming  results  obtained 
by  one  Fourier  method  into  its  equivalent  representation  under  another  form. 
We  thus  permit  conversions  (without  disturbing  questions  about  suitability 
and  applicability)  between  different  forms  of  these  Fourier  methods  used 
in  different  parts  of  the  literature.  In  particular,  the  turbulence  literature 
uses  significantly  different  nomenclature  and  transform  techniques  from  the 
propagation  literature.  Suitable  transforms  arc  possible  when  the  turbulence 
results  have  been  translated  into  propagation  literature  form.  These  processes 
constitute  the  bulk  of  chapters  2  and  3.  In  chapter  4  these  techniques  are  applied 
to  the  problems  of  phase  and  deflector  screen  propagation. 

Chapter  4  contains  an  overview  of  the  phase  screen  generation  and  derivation 
process.  Phase  screens  have  been  used  in  the  past  largely  in  beam-wave 
applications.  Our  interest  in  image  modification  requires  working  with  large 
area  effects.  To  support  this  goal,  a  less  computationally  intensive  deflector 
screen  method  is  introduced.  This  method  extends  prior  refractive  raytracing 
and  mirage  imaging  techniques  that  handled  only  vertical  inhomogeneities  in 
the  refractive  index  structure  for  fully  three-dimensional  effects. 


Conclusions 

In  this  report,  the  mathematics  of  several  key  concepts  are  developed  in  support 
of  simulating  turbulence  effects  on  image  propagation.  As  a  result  of  these 
developments  the  underlying  means  of  translating  between  different  forms  of  the 
Fourier  transform  was  shown.  This  development  permitted  the  writing  of  the 
turbulence  spectrum  in  terms  that  can  be  used  in  computer  representation  of 
turbulence  effects.  In  combination  with  an  analysis  of  both  the  phase  screen 
methodology  and  the  deflector  screen  approach,  we  showed  that  large  scale 
turbulence  can  have  significant  effects  on  image  distortion.  A  means  was  also 
discovered  for  handling  finite  sized  incoherent  source  regions  emitting  energy 
that  is  received  at  a  circular  aperture.  The  resulting  analysis  showed  that  high 
frequency  turbulence  has  little  effect  on  image  shape  distortion.  The  result  of 
these  analyses  is  a  method  for  distorting  imagery  for  heat  boil  type  turbulence 
distortion  effects.  These  effects  essentially  consist  of  a  combination  of  angle- 
of-arrival  fluctuations  which  distort  apparent  object  location  and  isoplanatism 
which  causes  different  portions  of  the  scene  to  distort  in  different  directions. 
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1.  Introduction 


Recent  developments  in  optics  have  resulted  in  sensors  of  increasing  sensitivity 
at  infrared  (TR)  wavelengths  with  significantly  reduced  noise  that  often  limited 
earlier  devices.  This  has  caused  an  increased  awareness  within  the  TR,  sensor 
development  community  that  turbulence  effects  should  be  included  in  sensor 
performance  models. 

Tn  particular,  turbulent  image  distortion  can  be  significant  for  TR,  wavelength 
imagery.  However,  the  modeling  tools  to  simulate  distortion  effects  have 
been  nonexistent.  One  reason  for  this  deficiency  is  that  prior  research  in  the 
turbulence  area  has  focused  on  astronomical  problems,  involving  plane  wave 
propagation,  narrow  fields  of  view,  large  aperture  sizes,  and  low  turbulence.  By 
comparison,  the  Army  problem  of  point-to-point  observation  in  which  observers 
attempt  to  acquire  and  track  targets  near  the  ground  usually  involves  spherical 
wave  propagation,  wide  fields  of  view,  aperture  sizes  close  to  the  Fresnel  zone  in 
diameter,  and  high  turbulence.  For  point-to-point  observations  with  both  object 
and  observer  located  near  the  earth’s  surface,  significant  optical  turbulence 
effects  may  be  encountered.  For  astronomical  problems,  the  observation  site  is 
usually  selected  to  avoid  as  much  turbulence  as  possible.  Often  the  observatory 
is  thermally  controlled  to  avoid  temperature  differences  between  the  outside 
atmosphere  and  the  entrance  pupil  to  the  telescope.  Such  luxuries  are  not 
possible  for  Army  operations. 

Flat  desert  conditions,  for  example,  can  cause  extreme  turbulence  conditions. 
When  discussing  turbulence  viewing  with  Army  and  Marine  officers  at  Ft.  Knox, 
KY,  some  years  ago,  a  comment  was  made  that  targets  more  distant  than 
approximately  1  km  were  simply  not  engaged  because  the  amount  of  turbulent 
distortion  and  blur  present  made  it  nearly  impossible  to  identify  targets  at  longer 
ranges.  At  IR.  wavelengths,  turbulence  effects  are  considerably  less  than  at 
visible  wavelengths  and  so  there  is  the  hope  that  IR  sensors  will  ultimately 
provide  better  performance  under  high  turbulence  conditions.  The  question  is 
how  much  better?  And  can  certain  tradeoffs  or  adaptive  optics  further  improve 
sensor  performance?  The  only  way  to  quantitatively  answer  such  questions  is 
through  accurate  modeling  of  turbulence  effects  and  human  perception  testing 
to  determine  how  observers  interpret  distorted  or  partially  corrected  signals. 
Modeling  synthetic  scenes  including  both  system  and  turbulence  effects  means 
the  ability  to  perform  perception  testing  under  controlled  conditions.  Planned 
adaptive  optics  augmentation  of  sensors  can  also  be  tested.  The  addition  of 
adaptive  optics  merely  makes  the  sensors  more  complicated  and  extends  the 
modeling  required  to  simulate  the  image  acquisition  process.  Such  a  system 
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would  still  have  to  be  evaluated  in  terms  of  its  ability  to  compensate  for 
turbulence  effects,  and  human  testing  would  still  be  required  to  determine  the 
efficacy  and  cost  effectiveness  of  such  methods. 

One  key  concern  in  determining  turbulence  effects  is  to  adequately  model 
the  differences  in  turbulence  effects  over  a  range  of  turbulence  strengths.  In 
this  report  we  discuss  methodologies  for  simulating  turbulence  within  image¬ 
rendering  software.  While  this  topic  is  not  new,  prior  efforts  have  been  primarily 
related  to  image  blur  effects  only.  Here,  the  focus  is  on  image  distortion. 

This  effort  has  strong  links  to  prior  efforts  to  understand  turbulence  effects  on 
beam  wave  propagation.  Unfortunately  that  literature  is  not  as  well  documented 
as  might  be  imagined.  While  several  studies  have  been  conducted  in  the  area  of 
beam  wave  propagation,  there  are  errors  within  the  standard  documentation. 
Later  papers  have  simply  repeated  these  errors,  while  others  have  provided 
correct  equations  without  providing  the  basis  for  these  results.  A  general 
presentation  of  turbulence  propagation  methods  is  thus  in  order,  a  presentation 
which  attempts  to  generate  seamless  and  consistent,  error-free  results,  while 
discussing  the  various  caveats  and  limitations  involved.  In  addition,  a  method 
based  on  raytracing  technology  is  discussed  that  should  provide  nearly  the  same 
fidelity  of  calculation  as  is  available  from  the  standard  phase  screen  methods 
without  the  significant  cost  in  computational  time  necessary  to  handle  the 
propagation  of  incoherent  light  at  the  level  of  tracking  the  propagation  of  field 
amplitudes. 

To  accomplish  this  analysis  task,  we  establish  the  nature  of  the  Fourier  transform 
and  its  relationship  to  other  Fourier  relations.  This  first  step  is  designed  to 
support  a  transformation  of  results  between  different,  but  related,  standards 
used  in  the  literature.  This  literature  crosses  several  disciplinary  boundaries, 
from  micrometeorology  to  statistics  to  optics  to  computational  methods.  Each  of 
these  disciplines  has  individualized  means  of  representing  Fourier  objects.  The 
separate  disciplines  must  be  examined  and  interrelated  to  smoothly  translate 
between  the  different  forms.  In  particular,  to  represent  the  Fourier  version 
of  the  refractive  index  power  spectrum  for  optical  turbulence,  one  must  be 
able  to  (1)  translate  between  two  forms  of  the  Fourier  transform,  and  (2) 
express  these  results  in  a  form  compatible  with  the  inverse  Fourier  transform 
using  a  Fast  Fourier  method.  This  process  involves  converting  the  Fourier 
transform  of  a  continuous  function  to  a  Fourier  series  coefficient  set  compatible 
with  the  periodic  functions  and  limited  frequency  regime  of  the  Fast  Fourier 
Transform.  This  form  must  be  compatible  with  a  propagation  code  that  handles 
sensor  optics  and  the  diffraction  propagation  between  turbulence  phase  screens. 
To  support  these  various  developments,  the  main  analyses  have  been  divided 
between  chapters:  Fourier  methods  in  chapter  2,  turbulence  representation  in 
chapter  3,  and  propagation  methods  in  chapter  4. 


2.  Fourier  Analysis 

In  this  chapter  we  discuss  Fourier  representations  of  various  types  of  functions. 
As  will  be  seen,  representing  these  functions  in  terms  of  sinusoidal  components 
yields  different  results  depending  on  the  nature  of  the  functions  themselves. 
These  differences  are  of  interest  because  of  the  various  ways  that  turbulence 
needs  to  be  or  has  been  represented  in  different  settings. 

In  particular,  in  previous  studies  of  optical  turbulence  structure  (cf.  Beland, 
1993)  the  refractive  index  spectrum  was  viewed  as  space  filling.  However,  we 
know  that  fluctuations  in  the  refractive  index  near  the  earth’s  surface  do  not 
extend  upward  significantly.  Further,  the  strength  of  turbulence  is  likely  a 
horizontally  varying  quantity.  Hence,  although  studies  of  turbulence  simulation 
discuss  the  nature  of  the  turbulence  in  terms  of  abstract,  infinite  fields  of 
turbulence,  when  treating  numerical  turbulence  within  finite  computers,  the 
turbulence  representation  is  always  finite.  Here  we  define  the  term  “field”  as 
a  multi-dimensionally  varying  random  function  which  obeys  certain  statistical 
properties.  Often,  in  describing  the  temporal  nature  of  turbulence  fields, 
Taylor’s  hypothesis  is  used  whereby  the  turbulence  is  considered  “frozen”  such 
that  it  advects  with  the  mean  wind  but  otherwise  does  not  change  its  spatial 
fluctuational  structure. 

In  terms  of  finite  computers,  the  turbulence  structure  is  specified,  limits  are 
applied,  and  fields  are  simulated  using  either  hardware  or  gridded  data  sets 
which  are  limited  in  extent.  Thus,  Fourier  series  representations  are  used  rather 
than  Fourier  transforms  (Andrews  and  Shivamoggi,  1999),  and  these  series  are 
represented  in  such  a  way  as  to  be  computed  by  fast  computational  methods. 

For  many  readers,  much  of  this  chapter  will  be  a  review,  but  the  information 
is  introduced  to  develop  a  common  framework  for  discussion  and  to  introduce 
terminology  and  the  particular  form  used  for  the  Fast  Fourier  Transform  (FFT). 

In  using  the  FFT  method,  the  modeled  volume  of  the  turbulent  region  is  assumed 
to  be  periodic  in  three  dimensions.  The  analysis,  then,  needs  to  consider  a 
three-dimensional  problem.  However,  our  analysis  will  begin  by  studying  a  one- 
dimensional  case  and  then  extending  the  results  to  three  dimensions.  The  x  axis 
is  chosen  for  the  analysis,  but  the  properties  discussed  will  be  the  same  on  each 
of  the  other  axes  following  changes  in  axis  labels  and  integration  limits. 
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2.1  Gaskill’s  Fourier  Transform 

We  begin  by  discussing  a  form  of  the  Fourier  transform  (Gaskill,  1978)  that  has 
symmetric  forward  and  inverse  transform  equations.  This  form  will  be  followed 
closely  when  considering  the  Fourier  series  representation  for  periodic  functions. 
The  results  of  these  considerations  will  then  be  used  in  relating  the  transforms 
currently  adopted  in  the  turbulence  literature  to  the  forms  needed  to  represent 
phase  screens  using  an  FFT  approach. 

2.1.1  Transform  Definition 

According  to  Gaskill  (1978),  the  Fourier  transform  is  defined  as 

OC 

F(0=  J  f(x)exp(-i2nx£)dx.  (1) 

—  OC 

One  can  then  recover  the  original  function,  /(#),  through  the  inverse  transform 

OC 

f(x)  =  J  F(f)  exp(*27ra:£)  d£.  (2) 

—  OO 

The  function  F(£)  is  called  the  Fourier  transform  of  the  function  f(x),  and 
similarly  the  function  f(x)  is  called  the  inverse  transform  of  F(£).  Note  the 
symmetry  involved  in  Gaskill’s  version  of  the  Fourier  transform.  Other  forms, 
including  those  used  in  much  of  the  turbulence  literature,  involve  external  2n 
factors  in  either  the  transform  or  inverse  transform  process.  (In  fact,  in  some 
of  the  literature,  the  transformed  functions  are  not  called  Fourier  transforms 
at  all.  Instead,  they  are  referred  to  as  characteristic  functions  (Papoulis,  1984; 
Panofsky  and  Dutton,  1984).) 

Under  GaskilPs  transform  definition,  a  dilated  and  shifted  function  transforms 
as 

/  H  exP (±*27rxoO  F(a£),  (3) 

where  =>  denotes  the  Fourier  transform  operation. 

2.1.2  Transforms  of  Selected  Functions 

Within  this  framework,  the  transform  properties  of  several  functions  are  of 
interest.  The  first  of  these  is  the  Gaussian  function: 

Gaus(x)  =  e~*x  .  (4) 

The  Gaussian  is  one  of  several  functions  which  are  invariant  under  GaskilPs 
transform:  T 

Gaus(or)  Gaus(£).  (5) 
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We  can  use  a  limiting  version  of  this  function  in  our  study  of  the  properties  of 
the  impulse  function,  also  commonly  known  as  the  Dirac  delta  function.*  The 
properties  of  the  Dirac  delta  function  flow  from  the  defining  equation 


I  f(y)S(y-x)dy  =  f(x)  (6) 

o 

for  suitably  smooth  functions  f(x).  This  definition  is  also  referred  to  as  the 
delta  function  sifting  property.  We  can  study  the  Fourier  transform  properties 
of  a  delta  function  using  a  limiting  form  of  the  Gaussian  function  in  combination 
with  the  dilation  properties  noted  in  equation  (3).  The  delta  function  can  be 
written  as 

<5(x)  =  lim  A  Gaus  .  (7) 

Ob — >0  | Cl |  V  (X  * 

The  area  under  the  function  on  the  right-hand  side  (RHS)  of  equation  (7)  is  unity 
for  all  nonzero  a.  Integration  of  this  function  times  f(x ),  as  in  equation  (6), 
yields  an  approximation  of  /( 0)  which  improves  as  u  — >  0.  In  the  limit,  this 
result  is  written  using  the  notation  in  equation  (6);  however,  strictly  speaking, 
the  limit  cannot  pass  through  the  integral  operation,  as  equation  (7)  does  not 
satisfy  the  condition  of  bounded  convergence  (see  the  previous  footnote).  One 
way  to  think  of  the  delta  function  in  practical  terms  is  as  a  sampling  function 
whose  operation  is  extremely  rapid  when  compared  to  the  rate  of  variations  of 
the  sampled  function  f(x). 

The  delta  function  also  features  a  scaling  property  such  that 

8  )  =  |a|  8(x  -  x0)-  (8) 

Finally,  based  on  the  dilation  rule  (equation  (3)),  we  can  transform  the  delta 
function  as 

5(x)  =$  lim  Gaus  (ax)  =  1.  (9) 

a— ►() 

The  <5  function  and  its  properties  are  necessary  to  understand  the  various  features 
and  properties  of  Gaskill’s  “comb”  function: 

(—  \  00 

- — —  j=|a|  8(x  —  xo  —  na).  (10) 

'  n=— oc 


*  A  rigorous  treatment  of  the  delta  function  involves  the  use  of  generalized 
functions  (Kolmogorov  and  Fomin,  1970).  The  approach  taken  here  views 
equation  (6)  as  a  notation  representing  the  linear  functional  Lx(f )  =  f(x)  rather 
than  an  actual  integral. 


The  name  for  this  function  (characteristic  of  other  Gaskill  function  names) 
evokes  its  visual  appearance.  In  this  case,  the  comb  consists  of  an  infinite  series 
of  delta  “tines”  occurring  at  separation  intervals  \a\  along  the  x  axis. 

The  comb  function  is  also  invariant  under  Fourier  transformation: 

comb(x)  comb(£).  (11) 

GaskilPs  (1978)  discussion  of  this  relation  is  rather  detailed  (pp.  205-206) 
and  will  not  be  repeated  here.  In  it,  Gaskill  admits  to  ignoring  “the 
Dirichlet  condition  that  prohibits  impulsive  behavior  in  the  function  to  be 
expanded.”  The  result  is  thus  not  a  rigorous  proof,  but  can  be  made  so 
using  careful  application  of  the  theory  of  generalized  functions  (Kolmogorov 
and  Fomin,  1970). 

Using  the  comb  function,  a  uniformly  spaced  sequence  of  delta  functions  will 
have  a  particularly  simple  form  in  frequency  space: 

OC  1 

7;  S(x  —  na)  =  pr  comb  (C'j  =*>comb(a£).  (12) 

n . =  —  rv^  '  * 


2.1.3  Convolution  and  Cross  Correlation  Operations 

Given  the  form  of  the  Fourier  transform  we  see  immediately  that  multiplication 
of  a  function  by  a  constant  scales  the  Fourier  transform  of  the  resulting  product 
by  the  same  constant.  We  have  already  considered  dilated  and  shifted  functions 
in  the  previous  section.  However,  there  are  two  other  linear  operations  of 
particular  interest  which  have  unique  Fourier  transform  properties.  These  are 
the  convolution  and  cross  correlation  operations. 

For  two  square  integrable  functions,  /(&*)  and  g(x)->  defined  over  the  x  axis,  we 
define  the  convolution  operation  as 


f(x)*g(x)  =  j  f(a)  g(x  —  a)  da.  (13) 

—  OG 

For  the  cross  correlation  operation  we  write 

OO  OC 

f(x)*g(x)=  J  f  (a  +  x)  g(a)  da  =  J  f  (a)  g(a  -  x)  da.  (14) 

— oc  — OC 

Gaskill  notes  (equation  (6.11))  that  while  f(x)  *  g(x)  =  g(x)  *  f(x),  the  cross 
correlation  operation  does  not  commute  (equations  (6.46)  and  (6.48)): 

f(x)  *  g(x)  =  f(x)  *  g(-x)  ^  g(x)  *  f(-x)  =  g(x)  *  f  (x).  (15) 
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Under  Fourier  transform  (Gaskill’s  Table  7.2)  these  operations  considerably 
simplify: 

f(x)  *g(x)  =,  F({j  G(-0- 

On  the  other  hand,  products  of  functions  transform  as 


2.1.4  Extension  of  Gaskill’s  Methods 

Gaskill  extends  his  one-dimensional  analysis  to  two  dimensions  in  order  to  use 
this  theory  in  Fourier  optics  applications.  We  can  extend  these  results  to  three 
dimensions  in  a  similar  manner.  For  example,  a  three  dimensional  version  of  the 
delta  function  becomes 

S(x ,  y,  z)  =  5{x)  5{y)  <5(z),  (18) 

Similar  extensions  exist  for  other  one-dimensional  functions. 

A  final  feature  of  note  is  a  consideration  of  the  units  by  which  the  original 
functions  are  represented  versus  the  units  of  the  transformed  functions.  As 
we  shall  see,  the  units  associated  with  series  coefficients  in  the  Fourier  series 
representation  of  a  function  are  the  same  units  as  the  original  function,  whereas 
the  units  of  the  transform  of  a  function  f(x )  consist  of  the  product  of  the 
original  units  of  /  (x)  times  the  units  associated  with  the  variable  x.  Integration 
with  respect  to  the  variable  £,  whose  units  are  the  inverse  of  x.  removes  this 
dependence  during  the  inverse  transform  process. 

2.2  The  Fourier  Series 

In  the  representation  of  periodic  functions,  p{x),  a  Fourier  series  is  usually 
employed.  Consider  a  real- valued  function,  p(x),  which  is  periodic  over  intervals 
of  length  X : 

p(x)  =  p(x  +  ml),  (19) 

where  m  is  any  integer. 


7 


2.2.1  Transform  Definition 

Let  the  Fourier  series  representation  of  this  function  be  given  by 


p{x)  =  qp  +  2y 


( li  COS 


1=1  L 


/  27 Txl  \ 

V"*“/ 


+  bi  sin 


(¥) 


(20) 


This  form  was  chosen  to  closely  track  Gaskill's  Fourier  method.  Based  on 
orthogonality  relations  among  trigonometric  functions,  we  can  evaluate  the 
coefficients  a*  and  bi  as 

X/2 

(  2i vxl  \ 


ai 


=  Y  l  ^)cos(^) 


dx; 


-X/2 

X/2 


1  f  ,  .  .  /  2ttxI  \  , 

bl  =  x  /  p ^ sm  (~x~ )  dx' 


(21) 


(22) 


-X/2 


We  choose  the  region  of  integration  symmetrically  about  the  origin  for  reasons 
that  will  become  clear  later  in  the  discussion. 

The  factor  2  appearing  in  front  of  the  summation  in  equation  (20)  can  be 
removed  by  translating  these  results  into  complex  form  using  Euler’s  formula: 

el9  =  cos(0)  +  i  sin(0),  e~7^  =  cos(0)  -  %  sin(0),  (23) 

where  i  is  the  imaginary  root,  i  =  >/— I-  Rearranging,  we  have, 


cos (6)  = 


er 9  +  e  lQ 


jo 


-iO 


sin(0)  = 


(24) 


2  J  V7  2 i 

The  introduction  of  this  formula  into  equation  (20)  leads  to  the  complex  form 
of  the  Fourier  series,  with  the  terms 

2  [a,  cos  (27TX&)  +  bt  sin  (2wa*,)]  =  Pi  +  Pt*  e~i2^‘ ,  (25) 

where  =  l/X  and 

Pi  =  at  -ibp,  Pi  =  ai  +  ibi.  (26) 

A  superscripted  asterisk  denotes  the  complex  conjugate  operation. 

The  complex  form  of  the  Fourier  series  is  then 


p(x)  =  E  Pi  exp  (*27r^jx) , 

l=  —  OO 

X/2 

Pi  =  -t;  J  p(x)  exp  (— i2n£ix)  dx. 

-X/2 


(27) 


(28) 


Because  p{x)  is  a  real-valued  function,  the  coefficients  Pi  obey  the  Hcrmitian 
property  (P_j  =  P(*),  as  may  be  inferred  from  equations  (25)  and  (26)  above. 
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2.2.2  Fourier  Series  Convolution  and  Cross  Correlation 

Two  important  related  operations  involving  the  Fourier  series  are  again 
convolution  and  cross  correlation.  For  the  convolution  operation  we  write 

X/2 

c(x)  =  f(x)  *  g(x)  =  ~  j  f(x')g(x-x')dx'.  (29) 

-X/2 

The  cross  correlation  operation  is  given  by 

X/2 

q(x)  =  f(x)  *  g(x)  =  ^  J  f(x'  +  x )  g(x ')  dx' .  (30) 

- X/2 

This  function  becomes  an  autocorrelation  if  g(x)  is  replaced  by  f(x).  Notice 
that  these  operations  avoid  introducing  the  units  of  x  into  the  resulting  function, 
unlike  the  continuous  versions  in  equations  (13)  and  (14). 

We  are  interested  in  how  these  two  functions  transform:  Assume  f(x)  has  a 
Fourier  scries  with  coefficients  Fj,  and  g(x)  has  a  scries  with  coefficients  Gi.  We 
represent  the  Fourier  transform  of  the  cross  correlation  by  a  Fourier  series  with 
coefficients  Qc 


X/2 

Qi  =  -1  J  q(x)  exp  (— i27r  x  f  j)  dx 


X/2  T  X/2 

=T  j  JL  j  f(x'  +  x)g{x')dx'  exp  {-i2n  x  £i)  dx. 

-X/2  \_  -X/2 

Multiplying  by  1  =  (;xj)  { v 2 tt  x'  l/i)  exp  (— i2-n  x'  £j),  the  order  of  integration  can 
be  reversed,  and  a  new  variable  x"  =  x'  +  x  defined.  Holding  x'  constant  in 
the  outer  integral,  we  wish  to  integrate  over  x  in  the  inner  integral.  Since  x' 
is  constant,  dx  — ¥  dx"  and  we  may  integrate  the  inner  integral  in  terms  of  x" , 
leading  to 

x/2  if;',  X/2 

Qi  =  -1  I  g{x')  ^  I  f(x")  exp(—i2irx"£i)  dx"  exp  {i2irx'^i)  dx' .  (32) 

—X/2  [  -X/2 

The  inner  integral  may  then  be  evaluated  to  the  constant  FJ,  which  can  then  be 
factored  out  of  the  remaining  integral. 

X/2 

Qi  —  Fi  -1  J  g(x')  cxp[i2nx' £i)  dx'  =  Fx  G-i  =  Ft  G? .  (33) 


9 


~  -  2 

In  the  case  of  an  autocorrelation,  note  that  Ft  F*  —  Ft  is  nonnegative. 

The  procedure  for  evaluating  convolution  series  coefficients,  C*.  is  similar  and 
results  in  series  elements 

Ct  =  FtGh  (34) 

These  forms  are  relatively  close  to  GaskilPs  results  given  in  equation  (16). 
There,  the  transform  of  the  convolution  equalled  the  product  of  the  individual 
transforms.  Here,  the  convolution  coefficients  consist  of  products  of  pairs  of 
appropriate  coefficients  of  the  individual  transforms. 

Let  us  now  consider  the  extension  of  these  results  to  three  dimensions.  Define 
the  form  of  a  transform  for  the  function  f{x ,  y,  z)  as 

oo  oo  oo 

f(s)  =  Y  Y  Y  hm,nexp(i2Trs-ii,m,n).  (35) 

l=— oc  m=— oc  n=— oc 

The  Fourier  series  coefficients  are  evaluated  via 

X/2  Y/  2  Z/ 2 

Fi, m,n=xyz  f  dx  f  dV  f  dzf(s)exp(-i2irs-Zl,m>n),  (36) 

-X/2  -Y/2  -Z/2 

where  the  two  vectors  are  given  by 

-  A  m  ^  n  f  /oryX 

s  =  xi+,yj  +  zk,  Ci,m,n  =  x1  +  Y 3  +  Zk’  37 


where  «,  and  k  are  unit  vectors  oriented  along  the  x ,  y,  and  z  axes,  respectively. 
The  dot  product  operation  inside  the  exponentials  is  thus  shorthand  for 


_  p  _lx  my 
s  *  —  v  *r  v 


Since  the  functions  we  will  be  considering  will  be  real- valued,  the  resulting  three- 
dimensional  array  of  points  will  exhibit  a  three-dimensional  Hermitian 


property: 


F—l^—rn,—n 


(39) 


Also,  when  the  signs  associated  with  the  subscripted  indices  are  not  specifically 
indicated,  it  will  be  expedient  to  write 


We  may  imagine  that  the  Fourier  series  coefficients  Fimn  are  associated^  with 
specific  points,  &mn,  in  a  frequency  space  denoted  by  general  vectors  £.  As 
the  size  of  the  periodic  volume  V  =  XYZ  expands,  volume  dimensions  X, 
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Y,  and  Z  expand,  resulting  in  closer  spacing  of  the  sample  points  within  the 
frequency  space.  There  are  obvious  advantages  to  closer  spacings  of  points  as 
this  enhances  the  resolution  of  the  transformed  function.  As  the  number  of 
sample  points  increases,  the  magnitude  of  each  sample  point  simultaneously 
decreases  according  to  the  factor  l/V  =  1/{XYZ).  The  decreasing  magnitude 
of  individual  sample  points  is  thus  exactly  compensated  by  the  increasing  density 
of  sample  points. 

2. 2. 3  Transforms  of  Periodic  Functions 

In  studying  the  connection  between  Gaskill’s  Fourier  transform  methodology  and 
the  Fourier  series  representations  just  discussed,  we  now  consider  how  Gaskill’s 
Fourier  methods  handle  periodic  functions.  To  begin  this  discussion,  consider 
the  Fourier  transform  of  a  smoothly  varying  function,  f(x),  convolved  with  a 
delta  function.  From  the  rules  derived  in  section  2.1  we  have 

5(x)*f(x)^lxF(0  =  F(0.  (41) 

Thus,  the  delta  function  acts  as  the  identity  operator  of  the  convolution 
operation,  just  as  zero  is  the  identity  operator  under  addition  and  one  is  the 
identity  operator  of  multiplication. 

Consider  then,  the  convolution  of  f(x)  with  an  offset  5  function: 

S(x-x0)  *  exp(— * 27tx0£)  F(£)^>  f(x-r,o).  (42) 

We  find,  then,  that  convolution  of  f(x)  with  the  offset  delta  function  is  equivalent 
to  positionally  shifting  the  original  function  f(x)  by  the  distance  rro-  For  suitable 
functions  f{x),  satisfying 

b 

j  |/(x)|  dx  <  oo,  \f(x)\  <  oo,  (43) 

a 

where  a  and  b  are  limits  such  that  f(x)  =  0  for  x  <  a  and  x  >  b,  we  then  find 
the  convolution  of  f  (x)  with  an  appropriately  scaled  version  of  Gaskill’s  comb 
function  yields  the  periodic  function.  p(x), 

p(x)=p(x  +  X)  =  f(x)*j^  comb  (44) 

For  the  sake  of  simplicity  let  us  assume  that  f(x )  =  0  for  all  |x|  >  X/2.*  Thus 
individual  copies  of  f(x)  spawned  by  the  convolution  with  the  comb  function 
will  not  spatially  overlap. 

*  Technically,  a  and  b  could  occupy  any  interval  along  the  real  axis  and  we  need 
not  have  b  —  a  <  X.  The  simplifications  involves  establishing  symmetric  limits 
for  a  and  b  about  the  origin. 
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Prom  equation  (16),  we  can  see  that  p(x)  Fourier  transforms  into 

f{x)  *  -i:  comb  (-^)  F(£)  comb(X£).  (45) 

Using  the  definition  of  the  comb  function,  we  may  then  write 

P{X)^P( &  =  (46) 

Note  that  while  F(£)  was  defined  at  all  frequencies,  the  function  P(£)  is  only 
nonzero  at  unique  frequencies,  separated  by  intervals  =  1/X.  Further  note 
that  because  f(x)  equals  zero  outside  the  bounded  region  —X/2  <  x  <  X/2,  we 
may  rewrite  the  integral  for  F(£ )  as 

oo  X/2 

jF(£)  =  J  f(x)exp(—i2nx£)dx  =  j  f(x)  exp(— •/'  2ir  x£)  dx.  (47) 

-oc  -X/2 


Comparison  between  this  result  and  equation  (28)  reveals  a  correspondence 
between  the  Fourier  transform  representation  of  a  periodic  function  using 
Gaskill’s  Fourier  transform,  and  the  Fourier  series  representation  described  in 
section  2.1.  The  two  transforms  are  related  through  the  equation 

F(l/X) 


X 


(48) 


This  can  be  shown  by  deriving  a  formula  for  p(x)  starting  with  Gaskill’s  inverse 
transform: 


CO 

p(x)  =  /  P(£)  exp(*27rx£)  d£ 


OO 

S  1c)  exp(*27ra;^  ^ 


OC 


(49) 


^  F(l/X )  /  l 

-  22  — ^  -  exp  \  i2ttx 


l~  —  OC 
OC 


X 


X 


=  ^2  exp  (* 2;r  x  ■ 


/=— oc 


X 


The  final  form  of  this  equation  corresponds  to  the  complex  form  of  the  Fourier 
series  (equation  (27))  using  the  translation  law  contained  in  equation  (48)  above. 
We  thus  have  a  direct  method  of  comparison  and  inter-transformation  between 
two  different  Fourier  methods  for  periodic  functions. 

This  is  a  key  result-the  Fourier  series  is  merely  a  representation  of  the  Fourier 
transform  when  considering  periodic  functions. 


2.3  The  Fast  Fourier  Transform 


The  FFT  is  simply  a  method  of  calculating  a  discrete  Fourier  transform  that 
uses  a  series  of  acceleration  techniques  originally  developed  by  Danielson  and 
Runga  that  were  later  popularized  by  Cooley  and  Tukey  (cf.  Ludeman,  1986). 
Ludeman’s  treatment  of  this  technique  is  rather  thorough  and  is  followed  closely 
in  this  text.  Another  standard  treatment  is  found  in  Press  et  al.  (1992),  but 
involves  a  different  sign  in  the  complex  argument.  Regardless  of  the  treatment, 
the  key  feature  of  the  FFT  is  that  a  method  has  been  found  whereby  the  number 
of  complex  multiplications  required  to  transform  a  data  set  can  be  significantly 
reduced  by  processing  the  data  in  a  specific  series  of  stages.  The  amount  of  time 
savings  is  proportional  to  the  data  set  size:  The  FFT  process  requires  on  the 
order  of  N  In  N  complex  multiplications  compared  to  N2  such  multiplications 
for  a  straightforward  discrete  Fourier  transform.  As  a  result  of  this  reduction, 
nearly  all  transform  calculations  performed  digitally  utilize  the  FFT  procedure. 

The  reason  for  our  present  interest  lies  in  standard  descriptions  of  refractive 
turbulence,  which  are  always  in  the  form  of  a  power  spectrum  (a  Fourier 
frequency  domain  representation).  These  spectra  form  components  of  functions 
(described  in  sections  4.3  and  4.4),  which  must  be  inverse  transformed  using 
the  FFT.  Difficulties  arise  when  one  attempts  to  translate  the  power  spectrum 
into  a  digital  representation  compatible  with  the  FFT.  This  section  lays  the 
groundwork  for  that  translation  task. 

The  FFT  process  itself  is  a  fast  version  of  a  discrete  Fourier  transform.  And 
the  discrete  Fourier  transform  is  simply  a  band-limited  version  of  a  Fourier 
series.  The  reason  for  the  band  limit  is  due  to  the  finite  number  of  samples 
in  a  digital  data  set.  We  let  N  be  the  number  of  samples  of  data  in  our  data 
set.  These  samples  are  normally  assumed  to  be  point  measurements  separated 
from  one  another  by  some  interval  Ax.  Tn  our  case  we  are  interested  in  spatially 
separated  data  so  Ax  has  dimensions  of  distance.  Because  the  FFT  is  related  to 
the  Fourier  series  method,  it  is  assumed  that  the  data  is  periodic.  The  period 
used  is  X  =  N  Ax. 

For  compatibility  with  the  Fourier  transform  and  series  definitions,  let  us  assume 
that  the  N  samples  are  centered  around  an  index  value  of  zero,  corresponding 
to  the  point  x  —  0.  Numerically  these  points  are  assigned  indices  starting  at 
—N/2  +  1  and  continuing  through  index  N/2.  Call  the  index  variable  n,  and 
assume  that  each  of  the  N  samples  is  a  measurement  of  some  underlying  physical 
process,  f{x).  Let  us  call  the  nth  sample  of  this  process  /(n),  corresponding 
to  the  value  of  f(x)  at  location  x  =  n  Ax.  Assume  for  the  moment  that  this 
sample  is  taken  instantaneously. 
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We  should  then  be  able  to  approximate  the  computation  of  a  Fourier  series 
coefficient  (equation  (28)  via  the  following  process: 


X/2 

Fi  =  i  I  f{x)  exp  [  -*27t^x-]  dx 


-X/2 

1  N/ 2 

f(n)  exp[— i2irlA£nAx]  Ax 

n=— AT/2+1 
,  N/2 

=  x  E  /(")  exp 

n=-N/  2+1 
N/2 


,  1  X 


X 

X 


=  ^  E  exp 

n=— AT/2+1 

=  Hi)- 


,2nln 
'  X 


(50) 


Here,  the  integration  over  the  x  interval  has  been  approximated  using  a 
step  interval  Ax  =  X/N  and  assuming  that  the  sample  value  approximately 
represents  the  value  of  f(x)  over  each  interval. 


We  used  a  change  of  notation,  referring  to  the  approximations  using  the 
symbology  F(/),  because  the  sampled  versions  are  qualitatively  different  from 
their  continuous  cousins.  For  one,  the  last  form  of  the  computational  equation 
involves  the  quantity  2i vln/N  in  the  exponent.  The  computation  of  the  F(l) 
coefficients  has  thus  lost  all  connection  to  the  period  length  X  dependence. 
Hence  it  is  necessary  to  do  external  book  keeping  to  keep  track  of  the  relationship 
between  the  indexed  frequency  results  and  their  physical  meanings.  Further, 
because  the  data  are  no  longer  continuous,  unlike  the  Fourier  series  coefficients, 
the  series  of  F(l )  coefficients  only  have  N  unique  values.  We  can  see  this  by 
looking  at  frequency  coefficient  F(l  +  N ).  We  find  that  it  has  the  same  value  as 
coefficient  F(l).  This  is  because 


exp 


2tx  (N  +  l)  n 
N 


—  exp  [—i  (27r  n)]  exp  —i 


/  2i xln 


exp 


27 xln 


)]■ 


V  N 


)] 


(51) 


since  exp[— *27rn]  =  1  for  any  n.  One  may  therefore  attempt  to  calculate  as 
many  spectral  coefficients  as  desired,  but  only  the  first  N  values  will  be  unique. 
All  others  will  be  periodic  duplicates  of  the  first  N  results. 


We  may  simplify  our  writing  of  the  transform  operation  by  the  introduction  of 
the  so-called  “twiddle”  factors,  W$: 

N/2 

Hl)= „  E  f(n)WN,  1  =  -N/2  +  1,...,  X/2;  (52) 

n=-N/2+l 


where 


WN  =  exp  (-*  27T  /N) ,  W$  =  exp  (^-i  ■  (53) 

The  range  of  l  values  has  also  been  chosen  symmetrically  about  the  origin  l  =  0. 
We  now  show  that  the  inverse  FFT  process  Ls  given  by 

N/  2 

/(«)=  E  n  =  —N/2  + 1.  ...,N/2.  (54) 

l=-N/2+l 


This  relationship  can  be  ascertained  as  correct  by  substituting  the  definition  of 
F(l )  into  the  equation  for  f(n ): 


/(»)  =  £F(!)e*p(«^), 


=  £ 


££/(«•)  «p(-^)|»p(< 

rt '  '  '  > 


27ml 


N  )  ’ 


(55) 


/(n7)  v- 

=  £Vl>p 


.  27r  l  (n  —  n7) 

z - - - - 

N 


Due  to  the  orthogonality  between  different  elements  of  the  complex  inner 
summation. 


N/2 

E  exp 


i=— JV/2+1 


2nl(n  —  n7) 
N 


(56) 


where, 

5<m = { o;  J;::  <«> 

Here,  6jn  is  the  Kronecker  delta  function  with  l  and  m  integer  arguments.  In 
effect,  equations  (56)  and  (57)  are  the  mathematical  expressions  of  a  geometrical 
argument.  The  N  components  of  the  summation  in  equation  (56)  represent 
the  N  complex  roots  of  the  equation  xN  =  1.  These  roots  are  arranged 
symmetrically  about  the  origin  in  the  complex  plane.  When  summed  these 
roots  cancel  one  another  out  in  all  but  the  case  where  N  =  0. 


Using  the  Kronecker  delta  notation  we  can  rewrite  the  RHS  of  equation  (55)  as 


E/V)^'  =/(n),  Q.E.D.  (58) 

n' 


Note  that  5"  has  the  same  sifting  properties  with  respect  to  summation  that 
the  Dirac  delta  function  has  with  respect  to  integration. 
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There  are,  however,  problems  that  arise  with  respect  to  the  FFT.  These 
primarily  focus  around  the  limited  nature  of  the  data  set  in  terms  of  the  number 
of  samples,  IV,  and  the  finite  interval,  Ax,  between  samples.  Due  to  the  limited 
number  of  samples  and  their  assumed  periodicity,  using  raw  data  one  often 
encounters  jump  discontinuities  at  the  ends  of  the  data  sets.  To  limit  the 
amount  of  spurious  noise  that  would  be  introduced  to  the  spectra  if  this  data 
was  transformed  in  raw  form,  the  data  is  nearly  always  “windowed”  (Harris, 
1978;  Ludcman,  1986).  The  windowing  process  involves  multiplying  the  data 
set  by  a  function  which  tapers  to  zero  at  both  ends,  resulting  in  an  imposed 
periodicity. 

The  second  effect  relates  to  the  finite  intervals  Ax  and  X.  Because  of  the 
finite  number  of  samples,  the  minimum  resolvable  frequency  is  1/X  and  that 
maximum  unique  frequency  is  N/(2X)  =  1/(2  Ax),  which  is  also  referred  to 
as  the  Nyquist  frequency.  Usually  the  Nyquist  frequency  is  used  to  establish  a 
requirement  for  low-pass  filtering  of  the  input  analog  signal  to  remove  higher 
frequency  content.  Signals  arriving  with  energy  content  at  higher  frequencies 
than  the  Nyquist  frequency  will  be  aliased  into  lower  frequency  portions  of  the 
computed  spectrum. 

In  our  problem,  windowing  is  not  necessarily  a  concern  because  we  begin  in  the 
transformed  domain  and  inverse  transform  to  the  real  domain.  The  resulting 
real  function  is  therefore  guaranteed  to  be  periodic.  However,  we  must  ensure 
that  the  size  of  the  domain  is  sufficient  to  adequately  characterize  the  turbulence 
structure,  including  outer  scale  effects  at  low  frequencies  and  inner  scale  effects  at 
high  frequencies.  Within  the  surface  boundary  layer,  the  length  scales  associated 
with  turbulence  extend  from  about  an  order  of  magnitude  smaller  than  inner 
scales  of  several  millimeters  (Tofsted,  1991)  to  outer  scales  measuring  up  to 
dozens  of  meters  (Tofsted,  2000).  Means  of  handling  these  requirements  are 
discussed  in  chapter  4  when  considering  generation  of  phase/deflector  screens. 


3.  Optical  Turbulence  Structure 

In  chapter  2  we  considered  three  related  transform  methods  useful  for  describing 
reah  periodic  real,  and  sampled  periodic  real  functions,  respectively.  We 
discovered  that  the  Fourier  transform  can  be  related  to  the  Fourier  series  of 
a  periodic  function  through  use  of  Dirac  delta  functions.  We  also  discussed  the 
similarity  in  form  between  the  FFT  and  the  Fourier  series.  However,  we  must 
now  develop  a  connection  between  results  contained  in  the  optical  turbulence 
literature  and  deflector  propagation  screens  discussed  in  chapter  4. 

Relating  these  two  literature  sources  is  not  simple.  The  FFT  method  assumes 
the  function  being  modeled  is  periodic  over  some  fundamental  volume  V. 
Conversely,  the  optical  turbulence  literature  generally  derives  results  involving 
expectation  values  for  aperiodic  functions  defined  over  all  space.  Hence,  rigorous 
evaluation  of  expectation  values  would  require  an  averaging  operation  evaluated 
over  infinite  volume. 

To  avoid  confusion  between  the  methods  described  in  this  chapter  and  those 
discussed  in  the  previous  chapter,  we  will  use  position  vector  f  instead  of  s. 
The  primary  quantity  of  interest  in  this  analysis  is  the  variable  n,  the  refractive 
index  in  the  earth’s  atmosphere.  One  very  useful  spatial  structure  property  of 
n  is  the  spatial  covariance  function,  Fn,  given  by  (Goodman.  1985) 

Fn  (5,  r2)  =  ([n(ri)  -  (n(fi))]  [n(f2)  -  (n(r2))]\ 

(59) 

=  (n{fi)n(f2)^  -  ^n(ri)^n(r2)^, 

where  angle  brackets  (())  represent  expectation  operations,  and  vectors  fi  and 
r2  represent  two  positions  in  a  presumably  infinite  medium. 

Modeling  the  volume  as  infinite  in  size  usually  leads  a  homogeneity  assumption 
which  I  now  introduce.  However,  I  make  this  approximation  while  recognizing 
that  it  has  limitations  and  may  need  to  be  altered  in  the  future.  The  problem  is 
that  homogeneous  turbulence  is  not  physically  realizable  (Hinzc,  1987).  Hinzc 
then  proceeds  to  considered  the  degree  of  damage  caused  by  maintaining  this 
assumption.  Hinze’s  following  commentary  observes  that,  far  from  damaging, 
inhomogeneities  are  necessary  to  maintain  the  level  of  the  turbulence  itself.  Due 
to  dissipation  at  small  scales,  energy  must  be  introduced  at  the  large  scales  to 
maintain  the  energy  cascade  across  the  inertial  subrange.  When  this  energy 
is  unavailable  the  turbulence  level  decays.  Though  a  significant  consideration 
from  the  point  of  view  of  models  to  describe  the  turbulence  itself,  the  amount 
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of  energy  needed  to  maintain  the  cascade  is  not  large  compared  to  the  energy 
maintained  in  the  field  itself,  so  the  fractional  error  is  relatively  small. 

Granted  these  problems  with  the  homogeneous  assumption,  we  nevertheless 
employ  it  to  avoid  major  complications  in  the  mathematical  analysis.  The 
homogeneous  turbulence  condition  involves  assuming  that  the  mean  refractive 
index  is  constant  with  position  and  that  the  structure  of  the  turbulence  is 
also  constant  with  position.  The  first  of  these  conditions  can  be  stated 
mathematically  as.  n0  =  (n(fi))  =  (n(f2))  in  equation  (59). 

The  second  condition  states  that  Tn  is  dependent  only  on  the  vector  difference 

U(r„r2)  =  rn(r),  r  =  r2-n.  (60) 

A  further  consequence  of  the  homogeneity  assumption  follows  when  we  parse 
the  refractive  index  into  its  mean  and  fluctuating  components: 

n(x,y,z)  =  n0 +  ni(x,y,z).  (61) 

From  the  definition  of  the  covariance  in  equation  (59),  as  long  as  the  mean 
refractive  index  is  constant  with  position,  we  can  write 

Tn(f)  =  (ni(fi)ni(fi  +r)).  (62) 

The  covariance  function  F„  is  thus  an  autocorrelation  statistically  averaged 
across  all  possible  realizations  of  the  «i(f)  field. 

Because  of  the  homogeneity  property  of  the  ri\  field,  and  because  its  second 
order  statistics  are  assumed  uniform  with  position,  we  can  show  that 

X/2  Y/  2  Z/2  , 

rn(^)  =  xYz  \  iff  ni^1)ni^®1  +  ) 5  (63) 

'-X/2-Y/2-Z/2 

where  X,  Y,  and  Z  are  the  dimensions  of  a  finite  volume. 

Following  methods  of  stochastic  integration  discussed  by  Stark  and  Woods  (1986) 
and  categories  of  random  processes  such  as  discussed  by  Goodman  (1985),  as 
long  as  the  ni  field  is  strictly  stationary  up  to  fourth  order  joint  probabilities, 
we  may  pass  the  expectation  operators  through  the  integration  operators  and 
obtain 


.  X/2  Y/2  Z/2  . 

(  /  /  /  ) 

'  V  /O  _V  /O  _  7  /•> 


-X/2-Y/2-Z/2 


=  [  f  [  (ni^ni(si  + $))  d$i  =  j  I  ^Tn(s)dsi 
=  Tn(s)  f  J  j  dsx  =  XYZTn{s). 


(64) 
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The  XY Z  volumes  then  cancel  and  equation  (63)  is  shown  to  hold. 

Let  us  now  write  the  volume  integration  associated  with  equation  (63)  as 


A(s)  =  Iff  ni(si)  ni(si  +  s)  dsi. 


(65) 


Rewriting  the  original  equation,  we  have 


r„(«)  = 


<d@> 

XYZ' 


(66) 


If  it  were  possible  to  extend  X,  Y.  and  Z  without  limit,  we  note  that  A(s)  then 
appears  as  an  autocorrelation  function  and  we  could  Fourier  transform  A  such 
that 

aw  =4  n(|)  =  mi-i).  (67) 

Ni(£)  here  would  represent  the  Fourier  transform  of  the  ni(s)  fluctuation  field. 
Unit  analysis  reveals  that  K  has  units  of  volume  squared  (m6)  while  Ay  has  units 
of  m3. 

Unfortunately,  this  procedure  is  not  possible.  The  resulting  integral  definition 
for  A  would  diverge  for  any  non-zero  value  of  Tn(s)  such  that  the  ratio 
(A(s))  /  XY  Z  remained  non-zero.  Because  of  this  limitation,  the  most  direct 
method  of  simulating  the  effects  of  turbulence  (by  directly  simulating  the 
fluctuation  field  n\  via  its  Fourier  transform)  is  unavailable. 

3.1  The  Refractive  Index  Spectrum 

It  is  in  the  homogeneous  form  that  we  normally  find  the  relationship  between 
the  covariance  and  the  refractive  index  power  spectrum,  $n: 


$«  («) 


Q  OC 

i)  III 


(r)  e*™  dr. 


(68) 


This  relationship  has  the  form  of  a  Fourier  transform,  but  it  is  immediately 
obvious  that  it  is  not  in  the  same  form  as  equation  (1).  The  difference  is  in  the 
change  of  conventions  between  the  standard  used  in  the  turbulence  community 
and  the  standard  used  in  writing  the  FFT  and  Gaskill’s  form  for  the  Fourier 
transform.  This  discrepancy  was  a  major  reason  for  writing  this  report  in  the 
first  place.  To  resolve  these  differences,  we  will  develop  the  refractive  index 
spectrum  according  to  the  standard  methods  used  in  the  turbulence  literature 
first,  and  then  find  the  means  of  representing  fl>n  in  a  form  compatible  with 
GaskilPs  Fourier  transform  in  section  3.3. 

The  quantity  it  is  called  the  spatial  frequency  vector.  Its  magnitude,  the 
spatial  frequency  k.  is  given  in  units  of  radians  per  meter.  The  refractive  index 
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covariance  can  be  recovered  from  the  refractive  index  power  spectrum  through 
the  inverse  transform  relationship 

oo 

rn(r)  =  jjj  <Mk)  e-ir^dK.  (69) 

—  OO 

The  power  spectrum,  4>n,  is  also  often  referred  to  in  statistics  literature  as  the 
characteristic  function  of  Tn  (Papoulis.  1984;  Panofsky  and  Dutton,  1984). 

As  a  further  simplification,  it  is  often  assumed  that  the  turbulence  is  isotropic. 
In  this  case  only  the  magnitudes,  r  and  k,  of  vectors  f  and  k,  respectively,  are 
significant.  By  integration  over  the  angular  variables  in  equations  (68)  and  (69), 
these  three-dimensional  integrations  can  be  simplified  to  the  one-dimensional 
relations  (cf.  Belaud,  1993): 

3  00 

*.<«)-(£)  W» 

o 

oo 

rn(r)  =  J  $w(k)  dn.  (71) 

o 


In  Tofstcd  (2000)  the  properties  of  the  outcr-scalc-infiucnccd  refractive  index 
spectrum  were  described  using  the  formula 


*»(«)  =  &C 


Ax  (. CaK f  £q1/3  (1-Tlr)  (W  41/3 


oi  17/6  r  ol 

1  +  (Ai  «)  1  +  (Cb  k)  J 


23/6 


(72) 


where  C2  is  the  refractive  index  structure  parameter,  having  units  of  m  2/3;  the 
constant  A\  was  found  to  have  the  value  8.2;  and  Ca  and  £(,  were  found  to  be 
related  to  the  outer  scale  through 

Ca  =  2.0741  L0;  £b  =  2.4767  L0.  (73) 

f3n  is  a  constant  of  integration. 


/3n  =  —  ^C5/6).  w  °-033> 

Pn  36  ir3/2  T(2/3) 

where  the  V  functions  here  refer  to  Euler’s  gamma  function,  defined  by 


oc 

T(a:)  =  J  e-<  tx~x  dt,  x  >  0, 


(74) 


(75) 


such  that  T(k  +  1)  =  kT(k)  =  k\ 


This  spectral  form  provides  a  more  accurate  model  of  the  turbulence  spectrum 
in  the  vicinity  of  the  outer  scale.  Unlike  the  misnamed  von  Karman  spectrum, 
which  approaches  a  nonzero  constant  value  at  zero  frequency,  this  spectrum 
correctly  approaches  zero  at  k  =  0,  consistent  with  von  Karman’s  original  intent 
(Hinze,  1987;  von  Karman,  1948). 


3.1.1  Outer  Scale  Definition 

In  Tofsted  (2000)  I  developed  a  more  realistic  refractive  index  spectrum  based 
on  data  collected  during  the  Kansas  1968  experiment.  Based  on  that  analysis  I 
set  the  outer  scale,  L0)  based  on  the  point  at  which  the  knee  of  the  new  spectrum 
significantly  departs  from  the  standard  Kolmogorov  spectrum  (Kolmogorov, 
1962),  symbolized  by  $nK-  To  identify  this  point,  let  k®  =  1/L0  equal  the 
frequency  associated  with  the  outer  scale.  Then  define  the  outer  scale  as  that 
point  in  the  spectrum,  where 


4>n  (K0)  =  («e>)  '  (76) 

The  Kolmogorov  spectrum  is  given  by  the  form 

$„*(«)=  A*  C£k-11/3.  (77) 

This  form  assumes  both  an  infinite  outer  scale  length  and  an  inner  scale  length 
of  zero.  Hence,  the  Kolmogorov  spectrum  is  applicable  for  the  inertial  subrange 
portion  of  the  spectrum,  and  this  is  normally  how  it  is  cavcatcd,  but  in  practical 
applications  involving  integrations  over  the  entire  spectrum  this  caveat  is  usually 
ignored. 


3.1.2  Covariance  Evaluation 


We  can  explicitly  evaluate  the  covariance  function  Fri  (r)  for  the  outer-scale- 
influenced  spectrum  through  direct  integration.  From  the  form  given  for  $n(fc) 
in  equation  (72),  it  is  expedient  to  identify  two  functions 


G\  (ra)  =  4tt 
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sin  ( ra  na) 
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To  solve  for  functions  G i  and  G<i.  one  could  use  integration  by  parts.  The 
term  k  dn/ (1  +  k2)v  can  be  integrated  directly,  through  the  change  of  variables 
x  =  1  +  k2.  By  this  approach  the  exponent  of  the  k  term  in  the  numerator 
is  reduced.  The  resulting  terms  can  then  be  integrated  using  Gradshteyn  and 
Ryzhik’s  (1980)  forms  3.771.2  and  3.771.5.  Waterloo’s  Maple  6  software  was 
used  to  initially  evaluate  these  integrals.  These  results  were  then  simplified  by 
combining  I-  type  Bessel  functions  to  produce  if -type  Bessel  functions: 


7r3/2  r¥3 

Gl  (r<l)  =  21/3r(17/6) 


[3ifi/3(r0)  -  raK2/3(ra)] ; 


(80) 
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3/2  1/3 

°2  ^  =  1721/3  r(l 7/6)  [(3r6  +  45)  K^Mrb)  -  26  rb  K2/3{rb )]  .  (81) 

Here,  T(r)  is  the  gamma  function  of  argument  r,  and  Ka(r )  is  the  modified 
Bessel  function  of  the  second  kind  of  order  a  and  argument  r  (cf.  Kreyszig,  1972, 
pg.  134). 

The  covariance  function  can  then  be  written 


rn(r)  =  /3nC I 


^/3Gi(£)+(1^i)£‘/3G2(£). 


(82) 


3.2  The  Refractive  Index  Structure  Function 

The  behavior  of  the  covariance  function  in  equation  (82)  can  be  compared  with 
the  behavior  of  the  Kolmogorov  spectrum  through  the  refractive  index  structure 
function,  Dn(f),  defined  as 

Dn(r)  =  ([«($)  -n(s  +  r)]2)  .  (83) 

Application  of  the  isotropic  and  homogeneity  conditions  reveals  that  Dn  can  be 
written  as  a  function  of  |r|  only.  Then,  Dn(r )  and  Tn(r)  are  shown  to  be  related 
through 

Dn{r)  =  2  [rn(0)  —  rn(r)] .  (84) 

The  use  of  Dn{r)  rather  than  Tn(r)  permits  comparison  between  the  outer- 
scale-influenced  spectrum  and  the  Kolmogorov  spectrum.  Different  covariance 
functions.  rn(r),  cannot  be  compared  directly  because  the  Kolmogorov 
spectrum  yields  a  singularity  at  Tn( 0),  indicating  its  overall  variance  is  infinite. 
Nevertheless,  Tatarski  (1961)  and  Clifford  (1978)  showed  that  the  structure 
function  for  the  Kolmogorov  spectrum  can  be  evaluated  directly  from  without 
evaluating  Tn  first.  Following  this  approach,  one  finds  the  Kolmogorov  form  of 
the  structure  function  to  be 


DnK(r)  =  C2nr2'\  (85) 

Normally  this  result  is  caveated  such  that  only  separation  distances  lQ  <  r  <  L0 
are  considered  valid.  Similar  restrictions  are  applied  to  the  range  of  valid  n 
values  in  the  Kolmogorov  spectrum  itself.  For  example,  Tatarski  (1961)  assigns 
Ko  and  Km ,  referring  to  outer  and  inner  scale  limits  such  that  Ko  <  n  <  Km. 
These  spectral  limits  define  the  bounds  of  the  so-called  inertial  subrange. 
However,  as  mentioned  in  section  3.1.1,  the  entire  spectrum  (0  <  k  <  oc)  is 
often  integrated  when  evaluating  the  effects  of  the  Kolmogorov  spectrum  in 
various  propagation  calculations.  It  is  for  this  reason  that  when  considering 
the  Kolmogorov  spectrum  in  the  remainder  of  this  document,  I  have  chosen 
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to  impose  the  full  consequences  corresponding  to  use  of  the  entire  Kolmogorov 
spectrum  to  highlight  the  deleterious  effects. 

The  first  of  these  effects  can  be  seen  when  comparing  plots  of  structure 
functions  based  on  the  outer-scale-influenced  and  Kolmogorov  spectra,  Dn(r ) 
and  Dnic(r)  respectively,  in  figure  1.  For  small  separation  distances  r  <C  L0, 
Dn(r )  and  DnK{r )  are  approximately  equal.  As  r  exceeds  L0  the  outer 
scale  limited  structure  function  levels  off  to  a  finite  maximum  value  2rn(0). 
Using  the  complete  Kolmogorov  spectrum,  we  find  that  rn#(0)  =  oo.  Hence, 
the  Kolmogorov  structure  function  increases  without  limit  with  increasing 
separation  distance  r. 


Normalized  Separation  (r/L0) 


Figure  1.  Comparison  between  Kolmogorov  and  outer-scale-influenced  structure 
functions. 


As  illustrated  in  figure  2,  as  the  size  of  the  outer  scale  increases,  so  does 
the  asymptotic  limit  of  Dn(r),  providing  a  close  approximation  between  the 
Kolmogorov  Duk  function  and  the  outcr-scalc-influcnccd  Dn  function  for  longer 
separation  distances.  This  figure  also  indicates  that  the  rise  distance  required  for 
Dn  (r)  to  reach  some  fixed  percentage  of  its  asymptotic  value  is  also  longer  as  La 
increases.  To  adequately  model  correlation  effects  in  a  propagation  simulation, 
one  must  contain  nearly  all  this  variation  within  the  fundamental  volume  V 
of  the  modeled  region.  Since  the  rise  distance  increases  with  increasing  La,  it 
naturally  implies  that  large  La  requires  larger  modeled  volumes  to  adequately 
simulate  the  influences  of  the  outer  scale  in  such  propagation  situations. 
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Figure  2.  Comparison  between  different  outer-scale-influenced  structure  functions 
with  varying  outer  scale  lengths. 

3.3  Transforming  the  Covariance  Function 

We  are  now  in  a  position  to  compare  the  Fourier-like  characteristic  functions 
used  in  the  standard  turbulence  theory  with  the  Fourier  methods  discussed 
in  chapter  2.  Due  to  the  difference  in  form  between  the  standard  form  of 
presentation  of  the  spectrum  as  opposed  to  its  means  of  implementation  in  the 
standard  form  for  the  FFT,  it  is  necessary  to  make  a  transformation  between 
the  two  forms.  In  particular,  we  need  to  compare  the  transformation  pair 
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with  the  transformation  pair 


oc 

*n$  =  flf  ^  e~i2*Hds,  (90) 

—  OO 

oo 

rn  (S)  =  ffl  *«  (l)  ei2*'*d£,  (91) 

— OO 

which  we  would  use  if  basing  our  methods  on  the  techniques  used  by  Gaskill. 

Previously,  we  have  used  s  and  r  to  denote  spatial  variables  in  the  two  different 
systems.  What  if  these  variables  actually  denote  the  same  quantity?  How  then 
are  functions  $n  and  related  to  one  another? 

— # 

Comparing  the  arguments  in  the  exponentials,  we  must  have  k  =  —  27r£.  We 
find,  however,  that  the  sign  difference  between  these  two  variables  is  immaterial 
because  both  and  must  be  radially  symmetric.  Thus,  the  magnitudes  of 
k  and  £  must  be  related  as  k  —  Making  this  substitution  in  equation  (90), 
we  find  that  \&n  is  equal  to  the  triple  integral  in  equation  (86),  leading  to  a 
relation  between  and  $n: 

'M£)  =  (27r)3*n(27r£).  (92) 

In  this  form,  the  \I>n  values  are  magnified  by  the  8n3  factor  because  in  £  space 
the  ’f'n  function  decays  proportionally  more  rapidly  compared  to  the  4>n  K-space 
representation.  The  difference  is  in  the  units  of  cycles  per  meter  in  £  space  versus 
the  radians  per  meter  units  in  k  space.  From  a  practical  viewpoint  these  results 
allow  us  to  avoid  converting  from  to  F„  and  then  to  T,, .  Instead,  we  can 
rescale  directly  from  <3>„  to  \kn. 

Having  established  the  relationship  between  and  4>n,  we  may  also  rewrite  the 
spherically  symmetric  forms  of  the  integral  relations  between  4>„(k)  and  rn(r) 
in  terms  of  1if,n(0  and  rn(s): 
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2 1T'% 

(94) 

This  spherically  symmetric  transform  pair  are  identical  (except  for  a  change  of 
variables)  for  both  forward  and  inverse  transforms.  These  are  thus  the  three- 
dimensional  equivalent  of  the  Hankel  transform  results  given  by  Gaskill  (1978, 
pg.  320)  that  were  also  self-reciprocal  for  forward  and  inverse  transforms  for 
cylindrically  symmetric  functions. 
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4.  Deflector  and  Phase  Screens 


This  chapter  considers  the  use  of  the  functions  F„  and  in  generating  what  will 
be  termed  deflector  and  phase  screens.  These  objects  are  usable  in  two  different 
approaches  to  modeling  propagation  through  optical  turbulence:  (1)  Deflector 
screens  are  useful  in  describing  the  shape  distorting  effects  of  turbulence  through 
a  model  of  positionally  varying  tilt  effects;  (2)  Phase  screens  model  turbulence 
effects  by  distorting  a  phase  front  and  propagating  the  distorted  phase  front 
between  successive  screens  via  diffraction  limited  propagation  equations.  This 
chapter  describes  both  propagation  methods  in  detail,  though  the  focus  is  on 
the  generation  of  the  screens  and  deflector  screens  in  particular. 

As  in  chapters  2  and  3,  gaps  also  exist  in  the  phase  screen  literature,  which 
are  again  addressed  through  the  derivations.  For  example,  a  key  assumption 
underlying  the  transformation  from  Fourier  transform  to  Fourier  series  results 
is  typically  handled  through  reference  to  prior  literature.  But  these  references 
are  to  papers  rather  far  afield  of  the  current  subject  matter  (e.g.  Shinozuka  and 
Jan,  1972;  Borgman,  1969).  Due  to  significant  differences  in  nomenclature  and 
symbology,  comparison  is  difficult.  Hence,  a  unified  presentation  is  sought  to 
restate  results  that  are  not  available  in  standard  references  on  this  topic.  Also, 
the  derivation  for  deflector  screens  is  completely  new. 

It  should  be  noted  here  that  the  overarching  goal  of  this  research  is  to  model 
the  propagation  of  a  scene  (an  undisturbed  original  image)  tlirough  synthetic 
turbulence.  The  model  for  this  propagation  assumes  a  radiance  pattern  (the 
scene)  is  being  emitted  from  a  source  region  (the  object  plane)  as  incoherent 
radiation  from  a  very  large  number  of  point  sources.  If  no  turbulence  were 
present,  the  energy  emitted  from  the  object  plane  would  propagate  over  a  path 
of  length  Z  to  an  imaging  system  of  aperture  diameter  D  and  focal  length  /. 
The  imaging  system  model  consists  of  the  single  aperture,  an  objective  lens, 
and  an  image  plane  where  the  received  image  is  detected  through  some  sort  of 
quantized  array.  This  array  samples  the  image  into  regions  called  pixels  (picture 
elements)  of  single-axis  angular  subtense  (f)p  radians  each.  Projecting  backward 
through  the  system  lens  and  over  path  length  Z ,  each  pixel  in  the  image  plane 
corresponds  to  energy  emitted  from  a  region  of  size  Z^>p  by  Z<f>p  in  the  object 
plane.  Even  for  very  narrow  field  of  view  systems,  each  such  source  region  must 
consist  of  a  very  large  number  of  independent  emitters,  in  support  of  the  concept 
that  the  signal  arises  as  an  incoherent  radiation  field. 

Based  on  this  model,  the  remainder  of  this  chapter  describes  how  phase  screens 
and  deflector  screens  can  permit  the  incorporation  of  turbulence  effects  to  the 
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propagation  scenario  just  described  for  diffraction  limited  optics.  Section  4.1 
provides  an  overview  of  diffraction  and  deflection  propagation  techniques.  In 
section  4.2  an  analysis  is  performed  for  light  passing  through  a  slab  of  atmosphere 
of  thickness  A.  Perturbations  in  the  crossing  time  for  light  entering  the  slab 
at  different  transverse  positions  can  then  be  modeled  in  terms  of  crossing  time 
correlation  spectra  described  in  section  4.3.  These  spectra  are  applied  directly  in 
the  diffraction  propagation  approach.  The  deflection  method  relies  on  a  further 
determination  of  beam  deflection  effects,  as  discussed  in  section  4.4. 

4.1  Propagation  Methodologies 

As  pointed  out  at  the  close  of  chapter  3.  functions  rn  and  cannot  be  directly 
inverse  transformed  to  generate  a  randomized  volume-based  n\  function.  If 
this  were  so,  randomized  n\  fields  could  be  generated  via  an  inverse  transform 
process.  A  propagation  technique  could  then  numerically  propagate  through 
this  volume.  As  will  be  seen  in  sections  4.1.1  and  4.1.2,  such  propagation  can 
be  modeled  in  two  fashions,  the  former  based  on  geometric  optics  and  the  latter 
based  on  diffraction  propagation. 

The  geometric  optics  approach  represents  an  extension  of  the  refractive 
raytracing  methodology  (Tofsted,  1989b)  developed  to  study  refractive  path 
bending  effects  bn  tank  gunnery.  This  method  is  akin  to  image  modification 
techniques  later  developed  to  study  the  effects  of  mirages  (cf.  Lehn,  2000; 
Lehn  et  al.?  1994;  Lehn  and  Friesen,  1992;  Sozou  and  Loizou,  1994)  on  image 
propagation.  Both  methods  considered  only  vertically  stratified  atmospheres 
and  only  dealt  with  the  mean  refractive  index  structure.  The  development  of 
deflector  screens,  which  covers  sections  4,2  through  4.4,  extends  these  effects  to 
three-dimensionally  varying  refractive  index  structures. 

The  development  of  the  deflector  screen  approach  is  based  in  large  part  on 
the  foundation  of  the  phase  screen  method  originally  developed  to  study  beam 
wave  propagation  through  turbulence  (cf.  Martin  and  Flatte.  1988;  Davis  and 
Walter,  1994).  The  original  developers  of  the  phase  screen  method  seem  to  be 
Fleck  et  al.  (1975),  but  their  published  analysis  contained  errors.  These  errors 
later  appeared  unchanged  in  the  open  literature  version  of  their  report,  Fleck  et 
al.  (1976).  Regardless,  in  sections  4.2  and  4.3  a  reanalysis  of  their  results  is 
given  which  correctly  derives  their  equations. 

In  general  the  term  screen  implies  a  two-dimensional  object.  In  effect  the 
use  of  screens  means  that  turbulence  effects  along  the  optical  path  are  being 
concentrated  into  planes  (or  screens)  of  interaction.  The  term  screen  is  used 
simply  because  a  mesh  of  sample  values  is  used  instead  of  a  continuous  function 
of  position.  Both  phase  screens  and  deflector  screens  have  these  qualities  in 
common,  as  their  derivations  are  similar. 

In  this  analysis  we  shall  see  the  other  ways  in  which  the  deflector  and  phase 
screen  approaches  are  similar.  However  in  its  full  form,  the  phase  screen 
(diffraction  propagation)  approach  should  produce  the  same  qualitative  results 
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as  the  deflector  screen  (geometrical  optics)  approach  while  yielding  more 
accurate  results.  The  difference  is  computational  cost.  To  be  accurate,  the 
phase  screen  approach  may  require  orders  of  magnitude  more  computational 
time  than  the  deflector  screen  approach,  yet  provide  only  marginally  different 
results.  As  pointed  out  by  Pried  (1982).  the  majority  of  effects  can  be  modeled 
using  geometrical  optics  methods.  Though  potentially  prohibitive,  phase  screen 
method  results  can  be  used  in  canonical  cases  to  provide  a  check  on  the  accuracy 
and  circumscribe  the  region  of  applicability  of  the  deflector  screen  approach. 

4-1.1  Geometric  Optics  Approach 

The  geometric  optics  method  is  similar  to  the  refractive  raytracing  calculations 
made  by  the  EOSAEL  (Electro-Optical  Systems  Atmospheric  Effects  Library) 
module  REFRAC  (Tofsted,  1987,  1989a,  1989b).  This  module  was  originally 
developed  to  study  the  effects  on  tank  gun  accuracy  of  vertical  refractive  index 
gradients  caused  by  surface  heating  and  cooling.  Here,  though,  the  technique 
is  extended  to  include  gradients  of  refractive  index  in  both  the  vertical  and 
horizontal  directions. 

Let  us  begin  by  defining  the  propagating  radiation  in  terms  of  a  small  angle 
approximation.  Under  this  approximation  we  assign  the  main  direction  of 
propagation  as  lying  along  the  positive  2  axis.  Let  us  identify  a  unit  vector 
Cl  as  the  direction  of  propagation,  consisting  of  the  three  components: 
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Since  the  propagation  is  almost  entirely  in  the  forward  direction  we  may 
approximate  the  vector  by 


Cl  ^  {OLx.  OCy ,  l)  ,  lor^l,  (o^yl  1. 


(96) 


Since  the  z  component  of  this  unit  vector  is  always  approximately  unity,  we  may 
ignore  this  component  and  study  the  tilt  induced  on  the  two-dimensional  vector, 
a  =  ( ax ,  ay). 


As  is  described  in  Tofsted  (1989b),  the  unit  Poynting  vector  associated  with 
a  propagating  wavefront  is  just  Cl.  We  may  model  the  tilts  induced  on  this 
wavefront  as  it  passes  through  a  non-uniform  atmosphere  by  considering  the 
gradients  of  the  refractive  index  oriented  perpendicular  to  this  wave.  Since  the 
wave  is  propagating  in  primarily  the  z  direction,  these  gradients  arc  located 
in  the  (x,  y )  plane  and  affect  the  vector  a.  The  effects  of  these  tilts  can  be 
described  by  the  equation 
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(97) 
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where  s  is  a  path  parameter  measuring  propagation  distance  of  a  point  along 
the  transverse  plane  in  the  direction  of  Cl,  and  V_l  is  the  gradient  operator  in 
the  transverse  plane.  This  equation  represents  an  approximation  of  the  results 
of  Born  and  Wolf  (1964,  pg.  Ill)  under  the  small  angle  assumption.  Born  and 
Wolf  relate  their  equation  to  the  eikonal  equation,  the  fundamental  equation  of 
geometrical  optics. 

In  Tofsted  (1989b)  the  refractive  index  was  assumed  to  only  vary  as  a  function  of 
height,  and  this  function  was  assumed  to  be  attached  to  the  surface  underlying 
the  optical  path.  Thus,  the  beam  tilts  induced  were  limited  to  only  vertical  tilt 
effects.  These  effects,  however,  were  only  evaluated  for  a  range  of  elevation 
angles.  Lehn  et.  al.  (2000,  1994,  1992)  independently  developed  means  of 
modifying  imagery  for  these  same  effects  by  remapping  pixels  from  an  original 
scene  to  raytraced  positions  in  an  image  plane. 

The  means  of  performing  this  remapping  involves  the  reciprocity  principle. 
According  to  reciprocity  (van  de  Hulst,  1980),  one  can  either  propagate  in  the 
forward  direction,  from  some  point  (po)  on  the  object  plane  in  a  particular  initial 
direction  (do)  and  attempt  to  shoot  (cf.  Press  et  al.,  1992)  a  ray  through  the 
entrance  pupil  of  the  receiver  system,  or  one  may  begin  at  the  image  plane 
and  propagate  a  ray  backwards  through  the  volume.  Eventually  the  backward 
propagating  ray  will  intersect  the  object  plane  at  a  point  and  in  a  direction  that 
would  propagate  through  the  receiver  entrance  pupil. 

The  use  of  deflector  screens  involves  the  replacement  of  the  continuous  form 
in  equation  (97)  by  a  discretized  version.  Propagation  of  rays  between 
subsequent  screens  is  modeled  as  following  unperturbed  straight  lines  similar  to 
O’Shea’s  (1985)  equations  for  tracing  between  subsequent  thin  lenses.  Beginning 
at  a  point  pm  at  the  mth  deflector  screen,  we  propagate  a  distance  Zm  to  the 
next  screen  in  the  direction  am : 

Pm+l  =  Pm  T  &m 

Passage  through  the  deflector  screen  then  perturbs  the  vector  a  according  to 

C^m+l  =  4"  Actm (Pm+l) )  (*^) 


where  1/n  «  1/no,  since  n0  »  n\  in  equation  (61).  Second,  the  transverse 
gradient  operator  and  the  path  integration  have  been  interchanged  from  the 
original  form  involving  an  integration  of  Vj_n(s)  along  the  path.  Third,  since 
no  is  constant,  only  n*  appears  in  the  gradient. 

Wavefront  tilt  over  a  path  increment  of  length  A  is  effectively  concentrated  into 
an  effect  that  is  applied  at  a  central  point  A/2  into  the  slab.  This  means  that  the 
propagation  distance  to  the  first  slab  is  Z\  =  A/2  while  the  distances  between 
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all  other  screens  are  Zn  =  A.  Lastly,  there  is  a  propagation  step  to  travel  from 
the  last  slab  through  the  receiver  entrance  aperture  and  this  distance  is  also 
A/2.  Of  course,  when  using  reciprocity,  these  steps  are  actually  computed  in 
reverse  order. 

The  deflector  screen  consists  of  a  function  which  prescribes  the  values  of  tilts, 
Act,  as  a  function  of  position,  s,  of  ray  arrival  at  the  screen.  In  essence, 
O’Shea’s  one-dimensional  thin-lenses,  ray-tilt  equation  is  replaced  here  by  a 
two-dimensional  ray-tilt  equation  dictated  by  the  statistics  of  the  atmospheric 
turbulence. 


4.1.2  Diffraction  Approach 

The  diffraction  approach  attempts  to  capture  amplitude  and  phase  distortions  in 
addition  to  angle-of-arrival  effects.  As  a  consequence,  the  propagated  field,  when 
passed  through  an  optical  system  should  produce  both  scintillation  and  blurring 
of  the  propagated  signal.  Such  improvements  do  not  come  cheaply,  but  under 
certain  conditions  may  be  unavoidable  for  accurate  calculations.  However,  the 
computational  burden  may  be  prohibitive.  Martin  and  Flatte  (1988)  required  a 
Cray  XMP  supercomputer  to  evaluate  their  results  for  a  limited  size  beam. 

Here,  we  follow  the  discussion  contained  in  Martin  and  Flatte  (1988),  but  modify 
the  nomenclature  to  avoid  conflicts  with  the  symbology  elsewhere  in  this  text. 
We  begin  with  radiation  of  wavelength  A  and  wavenumber  k  —  2ir/  A  propagating 
(again)  in  primarily  the  positive  z  direction.  The  complex  amplitude  of  this  field 
is  given  by 

(f>  =  X  exj>(—ikz).  (101) 

The  scalar  amplitude  satisfies  the  paraxial  wave  equation 

2ik^  +  V|X  +  2fe2mX  =  0.  (102) 

oz 


To  propagate  through  turbulence  with  this  equation,  Martin  and  Flatte  divided 
the  propagation  into  two  separate  effects.  First,  the  amplitude  field  propagates 
a  distance  A  through  uniform  media  in  a  diffraction-limited  fashion.  In  this 
intervening  space  n\  =  0.  Procedurally,  this  involves  Fourier  transforming  the 
paraxial  wave  equation  over  the  transverse  axes  x  and  y.  This  produces  the 
equation 

dX(a,  z)  _  ,47r2a2x(^  ^  (103) 


8z 


2k 


which  has  the  solution 


X(<r,  z  +  A_)  =  X(<r,  z)  exp 


-i- 


.47t2<t2  A 

2  k 


(104) 


The  use  of  the  nomenclature  z  +  A_  anticipates  the  next  step  used  to  handle 
turbulence  effects:  the  phase  screen  application. 
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As  the  term  phase  screen  implies,  the  turbulence  effects  are  concentrated  into 
so-called  “screens,”  or  planes,  of  turbulence.  In  this  case  the  screen  is  located  at 
z+ A.  At  the  phase  screen,  a  portion  of  the  wavefront  passing  through  the  screen 
at  position  s  =  (x,  y)  in  the  transverse  plane  is  phase  advanced  by  an  amount 
0(s).  Since  the  effects  of  turbulence  have  been  compressed  into  a  screen  of  zero 
width,  there  is  no  need  to  apply  diffraction  effects  for  propagation  through  the 
screen.  Thus  the  propagating  field  after  passing  through  the  phase  screen  is 
given  by 

X(s,  z  +  A+)  =  X(s,  z  +  A_)  exp[«0(s)] .  (105) 

Combining  effects,  we  have  a  four  step  propagating  process: 

1.  Discrete  Fourier  transform  the  wave  amplitude  function  X  across  the  plane 
perpendicular  to  z,  resulting  in  the  function  X. 

2.  Perform  diffraction  limited  propagation  of  X  from  the  z  plane  to  the 
z  +  A_  plane. 

3.  Use  an  inverse  discrete  Fourier  transform  to  recover  X.  at  the  z  +  A_ 
plane. 

4.  Apply  phase  screen  effects  to  the  function  X  at  z  +  A_  to  pass  it  through 
the  screen  to  the  z  +  A+  plane. 

This  process  is  repeated  until  all  phase  screens  have  been  crossed  and  the  receiver 
aperture  plane  reached. 

4.2  Slab  Crossing  Time  Analysis 

Note  that  in  the  previous  two  subsections  the  quantities  6(s)  and  Aa(s)  were 
critical  in  evaluating  the  effects  of  the  turbulence  upon  the  arriving  wave  field 
or  front.  The  derivation  of  these  quantities  is  the  key  topic  of  the  remainder  of 
this  chapter.  To  evaluate  these  quantities  requires  analysis  of  beam  propagation 
through  a  turbulent  layer  of  thickness  A.  However,  there  are  certain  assumptions 
being  made:  First,  we  are  assuming  that  rays  passing  through  the  slab  of 
thickness  A  are  travelling  approximately  parallel  to  one  another.  This  simplifies 
the  propagation  integrations,  but  ignores  the  divergence  of  rays  originating  from 
a  particular  region  associated  with  a  single  pixel  that  could  pass  through  the 
entrance  pupil  of  the  receiver.  But  the  approximation  is  nearly  correct  and  will 
be  used  as  a  basic  assumption  throughout  this  section. 

To  assess  9(s)  and  A 5(5)  we  need  to  first  consider  the  fluctuations  of  certain 
propagated  quantities  and  their  correlations  at  different  points  of  entry  through 
the  slab  of  thickness  A.  The  most  basic  of  these  quantities  is  sometimes  called 
the  optical  depth  (T  =  JA  n  dz)  but  we  prefer  to  analyze  the  time  required  for  a 

ray  to  cross  the  slab  (r  =  J0A(n/c)  dz).  Here,  c  is  the  speed  of  light  in  vacuum. 
n/c  is  thus  the  time  taken  by  light  to  move  a  unit  distance  in  air  of  refractive 
index  n . 
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The  analysis  of  beam  propagation  statistics  begins  by  considering  a  wavefront 
of  radiation.  This  same  construct  will  be  used  for  both  the  phase  analysis  and 
for  the  beam  tilt  analysis.  In  the  former  case  we  may  envision  a  propagating 
complex  amplitude  wavefront  X.  Tn  the  latter  case  we  consider  the  light  divided 
into  a  series  of  rays  (small  pencils  of  light,  each  having  a  small  overall  divergence 
in  solid  angle  and  a  small  transverse  extent  as  it  passes  through  a  turbulent  plane 
of  interest). 

Our  main  object  of  consideration  is  the  time  required  for  an  arbitrary  wavefront 
or  ray,  at  a  particular  transverse  position  s,  to  pass  through  a  slab  of  intervening 
atmosphere  of  thickness  A.  This  crossing  time  is  written  as 


r  = 

z 

The  starting  position  of  the  beam  is  (s,  z)  =  (x,  y,  z).  The  ending  position  is 
(s,  z  -j  A),  where  we  ignor  transverse  displacements. 

For  our  purposes  we  are  not  directly  interested  in  r.  Rather,  we  are  interested 
in  the  quantity 

Z  +  A 

St(s)=  I  ni^’  ^  dz',  (107) 

Z 

which  measures  the  perturbation  of  r  about  a  mean  transit  time.  Both  r  and 
6t  evaluations  assume  propagation  along  the  positive  z  axis.  Of  course,  ir  <  r 
since  n\  <C  no-  Nevertheless,  the  slight  deviations  of  n\  are  critical  in  assessing 
turbulence  effects. 

Paralleling  phase  screen  analyses  (cf.  Fleck  et  al.,  1976),  we  wish  to  evaluate;  a 
cross  correlation  function  between  the  time  of  transit  of  two  parallel  rays  crossing 
the  slab  at  points  S\  and  s2,  written  as 


(106) 


7"<St(s1)  S2)  = 


<5t(si)  6t(s2) 

z-{-A  z-\-A 


n i(si,  zi)  ni(s2,  z2) 


t  Z~t~Ljk  Z~v  i 

-  w 

'  Z  Z 

Z+ A  z+ A 

=  J  dzl  I  ^ 5!(^3) 


(108) 


Here,  the  expectation  operator  being  applied  to  the  entire  integral  is  passed 
through  the  integration  operators  and  applied  directly  to  the  ni  product, 
which  is  permissable  because  the  bounds  of  the  integration  are  unaffected 
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by  the  expectation  operator.  Under  a  homogeneous  and  isotropic  turbulence 
assumption,  we  may  use  Tn: 


m(*2,*2)^  =  i.  Tn  (7|5W2|2  +  (Z!-Z2)2)  -  (109) 

We  can  then  replace  Tn  with  the  inverse  transform  of  (equation  (91)).  Also, 
the  vector  arguments  of  r^r  can  be  replaced  by  the  difference  vector,  s  =  ( x  \  — 
X2,  yi  -  2/2): 


z-\-  A  z+ A  ex 


rsT(s)  =  ^  /  dz\  J  dz2  J j expfi  2ix  a  ■  s\  do 


Z  Z 

oc 


(110) 


I  exp[2  27T^(2l  -  Z2) ]  ’j'n  (VT^F+ff) 


where  a  represents  the  two-dimensional  vector  components  of  £  in  the  transverse 
plane,  and  £ z  represents  the  component  along  the  propagation  axis. 

This  form  permits  a  change  in  the  order  of  integration  such  that  the  z 
integrations  are  accomplished  first.  To  perform  these  integrations  we  will  revert 
to  the  variable  k2  —  2i r£z  to  make  the  mathematics  somewhat  simpler: 


z+ A 


z+ A 


J  dz\  J  dz2  exp[inz(zi  -  z2))  =  2 


1  ~  cos(kz  A) 


Kt 


(111) 


assuming  A  >  0.  Fleck  et  al.  (1976)  obtain  this  same  result  through  a  different 
method:  In  their  discussion  they  introduce  absolute  value  signs  around  the 
difference  Z\  —  z2.  Possibly  this  is  because  Tn  is  a  function  of  distance,  r,  which 
is  positive  definite,  but  this  should  not  matter  since  Fn  is  spherically  symmetric. 
Unfortunately,  the  introduction  of  the  absolute  value  complicates  the  integration 
such  that  they  generate  an  imaginary  term.  They  must  then  argue  away  this 
imaginary  term  before  they  can  arrive  at  a  result  equivalent  to  equation  (111). 

Fleck  et  al.  (1976)  then  assume  that  as  long  as  A  is  large  compared  to  correlation 
lengths  of  the  refractive  index  power  spectrum,  and  as  long  as  the  spectrum  is 
slowly  varying,  the  £2  integral  may  be  approximated  as 


(X 

/  (<?,  6)  —  [1  -  cos(/czA)]  diz 

•/  Kz 

-oo 

oo 

=  /*"(*  £)  4  u a)] 


d,Kz 
27 r 


(112) 


jMg,  0) 

27T 


oo 

/  —  [1  -  cos(kzA)]  dnz. 

J  Kz 
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Upon  a  change  of  variables,  u  —  kz  A.  Gradshteyn  and  Ryzhik’s  (1980)  result 
3.782  allows  us  to  solve  this  integral  directly,  as 

J  4  [!  -  ™s(k.zA)]  d,Kz  =  A  J  2[1~^°s(,x)]  du  =  2ttA.  (113) 

J  KZ  J  U 

—  OO  —OC 


The  approximation  used  in  equation  (112)  assumed  that  A  S>  L0.  We 
may  establish  some  sort  of  lower  limit  on  the  size  of  A  by  considering  the 
correlations  present  in  the  function  Tn.  Figure  3  shows  the  structure  of  the 
covariance  function  with  distance  normalized  relative  to  outer  scale  L0 .  As 
shown,  significant  correlations  extend  to  approximately  10  outer  scale  lengths. 
Beyond  10  outer  scale  lengths,  Tn  becomes  negative  but  remains  very  small. 
Since  these  correlations  extend  in  all  directions  about  any  point,  it  would  appear 
the  size  of  the  modeled  volume  should  be  at  least  20  L0  along  each  axis  direction. 


Figure  3.  Covariance  function  normalized  relative  to  its  amplitude  at  zero  and 
normalized  abscissa  measured  relative  to  the  outer  scale  length. 


The  size  of  the  modeled  volume  is  also  important  in  that  it  determines  the  lowest 
frequency  that  can  be  represented.  A  volume  that  is  not  large  enough  will  not 
be  able  to  capture  the  form  of  the  spectrum  at  sufficient  resolution.  We  will 
return  to  this  point  when  considering  the  resulting  spectrum  at  the  end  of  this 
analysis. 
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Finally,  the  evaluation  of  rgT  thus  far  has  reduced  the  problem  from  a  five-fold 
integration  to  a  two-dimensional  inverse  Fourier  transform: 


oc 

rsT(s)  =  ^2  U vt 'n{v,  0)  exp[«  2i TO  ■  s\  da 


//F#- 


0) 


(114) 


exp[i27T(j  *  s\  da. 


4.3  Crossing  Time  Fluctuation  Spectra 

Continuing  to  follow  the  arguments  of  Fleck  et  al.  (1975)  we  recognize 
that  A\£n(a,  0)/c2  is  the  Fourier  transform  of  rjr.  Let  us  then  consider 
equation  (108):  We  recognize  that  rsT  is  equivalent  to  the  expectation  of  an 
autocorrelation  operation,  where  the  autocorrelation  operation  occurs  over  space 
and  the  expectation  is  taken  over  all  possible  atmospheres  exhibiting  the  same 
turbulence  spectral  properties. 

The  random  field  Sr(s)  thus  defines  the  fluctuations  in  crossing  times  of  waves 
propagating  through  the  turbulent  slab  of  thickness  A.  Let  us  assume  that 
this  random  field  has  a  Fourier  transform  T(a).  The  transform  of  cross 
correlation  rsT  must  therefore  be  somehow  connected  to  the  average  of  the 
product  T(a)  T*(a )  since  this  product  represents  the  autocorrelation  of  the  field 
St(s ).  However,  progress  in  developing  this  relationship  is  only  possible  if  we 
assume  that  Sr  is  periodic. 

To  begin,  equation  (108)  can  be  written  in  terms  of  the  difference  vector  s :  If 
we  assume  that  6r  is  periodic  over  a  region  measuring  X  and  Y  in  size  in  the  s 
plane,  then  we  can  define  the  expectation  in  terms  of  a  spatial  average  over  this 
fundamental  periodic  region  as 


X/2  Y/2  . 

Sr($ 2  +  s)  Sr(s2)  ds2  \  .  (115) 

'  -X/2  -Y/2  ' 

In  this  form,  we  see  that  rsT  includes  an  autocorrelation  function  as  defined  for 
the  Fourier  series  form,  in  equation  (30).  Let  us  denote  by  R[m  the  Fourier 
series  coefficients  associated  with  the  periodic  version  of  r^T .  If  we  then  operate 
on  equation  (115)  to  generate  the  Fourier  series  coefficients  associated  with  this 
equation,  we  obtain 

Rim  =  (rhnfl*my  (lie) 


36 


But  we  know  that  taking  the  Fourier  transform  of  equation  (115)  results  in 


/ _  t>( „  _  \  _  A'I,n(<Xx,  Gyi  0) 

T 6t\S)  ^  ^y)  —  ^2 


(117) 


And,  from  extension  of  equation  (48)  to  two  dimensions,  we  have  the  following 
relationship  between  Rim  and  R(o): 


Rim  — 


Rib?) 


XY 


(118) 


We  now  have  Rim,  a  Fourier  scries  representation  of  rjT .  We  recall  from  section 
2.2  that  the  Fourier  series  coefficients  retain  the  same  units  as  the  transformed 
function.  Therefore,  Rim  has  units  of  s2.  The  standard  method  (see  Fleck  et 
al.,  1975;  Martin  and  Flatte,  1988;  Davis,  1994;  and  Yan  et  ah,  2000)  of 
generating  the  random  field  Sr  involves  generating  random  Tjm  values  that 
satisfy  equation  (116),  and  then  inverse  transforming.  Generation  of  random 
Tim  values  involves  coupling  an  envelope  function,  which  defines  the  frequency 
weighting,  with  a  white  noise  Gaussian  random  process.  The  resulting  form  for 
the  Tim,  coefficients  is  given  by 


Tim  —  Tim  Glm.  (119) 

Because  Tim  is  constant,  it  factors  out  of  the  expectation  operator  on  the  RHS 
of  equation  (116).  Because  of  this  property,  it  is  standard  to  define  G/m  such 
that 

(GimG*im)  =  l.  (120) 

Thus, 

Tim  =  \fhm  =  (121) 

The  Gim  coefficients  are  normally  represented  as  the  sum  of  two  zero-mean  unit- 
variance  Gaussian  random  variables  such  that  G;m  exhibit  Hermitian  properties: 

Gim  =  +  *a2(f,m)  (122) 

"  V2 

Fleck  et  ah  (1975)  identify  the  following  results  of  the  Hermitian  conditioning 
on  coefficients  ai(l,m)  and  ai (l,m): 

(di(Z,  m)  a.2(l,  m))  =  0; 

(ai(l,m) «!(/',  m')>  =  Sf  S%'  +  Sf  S~m' ;  (123) 

(0,2(1,  m)  02(1', m'))  ==  Sf  -  Sf  S~m' . 
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Hence,  we  have  random  variables  which  are  uncorrelated  from  point  to  point, 
and  which  obey  the  conditions 


ai(/,  m)  =  «i(— l,  —  m);  (124) 


a2(/,  rn)  —  —  a2(— l,  —m).  (125) 

As  Davis  (1994)  explains,  Gim  “represents  the  Fourier  transform  of  a  grid 
of  uncorrelated  Gaussian  distributed  random  numbers  representing  phases. 
The  proper  spatial  structure  function  corresponding  to  turbulence  statistics  is 
imposed  upon  the  random  phases. ..by  applying  a  filter....”  The  filter  function 
in  this  case  is  T)m. 

With  this  Hermitian  conditioning  we  can  see  that  G*lm  =  Hence  we  can 

directly  show  that  equation  (120)  is  valid: 


(Glm  G*lm)  —  (Glm  G-i-m) 

[di(Z,  m)  +  ia2(/,  m)]  [ai(— 1,  — m)  +  m2(— m)] 

Vf  y/2 

di(/,  m)ai(— /,  —  m)  \  j  iai(/,  m)a2(— /,  — m) 

2  /  +  \  2 
/*a2(Z,  m)ai(—l,  —m)  \  / 1 i2Q2(Z,m)a2(-/, -m) 

\  2  /+\  2 
=  i  +  i(o)  +  ?;(o)  +  ?:2(-i)  = 

2 


(126) 


These  results  can  then  be  combined  and  used  in  either  the  phase  screen 
or  deflector  screen  approaches.  Equation  (119)  can  be  evaluated  using 
equations  (121)  and  (122).  The  prior  discussion  in  the  section  then  reveals 


that 


(127) 


over  the  area  AY.  For  the  phase  screen  problem,  we  recognize  that  St  represents 
the  extra  time  it  takes  the  wavefront  at  position  s  to  propagate  across  distance 
A  compared  to  a  mean  propagation  time  r  =  A  «o/c.  As  a  result,  in  the  z  +  A 
plane,  after  the  passage  of  time  r,  at  position  s ,  instead  of  advancing  in  phase 
by  the  average  amount  k.A,  the  phase  only  advances  by  k,A  —  kcSr /no.  The 
phase  delay  is  thus 

0(.s)  —  k  —  St(s).  (128) 

n0 

This  phase  delay  appears  directly  in  the  screen  propagation  calculation  in 
equation  (105).  To  describe  the  equivalent  effect  as  measured  by  the  deflector 
screen  approach,  a  more  complicated  analysis  is  required.  This  analysis  is 
described  in  section  4.4. 


4.4  Beam  Deflection  Analysis 

We  now  consider  the  evaluation  of  raytracing  deflector  screens.  Prom  section  4.3 
we  know  how  to  evaluate  a  random  field  of  temporal  delays  5t.  We  now  extend 
that  analysis  to  determine  how  these  temporal  delays  translate  into  beam  tilts. 
Consider  figure  4  where  we  view  two  transversally  separated  points  associated 
with  two  points  along  a  wavefront  that  is  propagating  through  an  atmospheric 
slab  of  thickness  A.  Due  to  the  slight  differences  in  the  refractive  index  for 
the  two  paths,  the  transit  time  for  each  path  will  be  slightly  different.  For 
argument’s  sake  let  us  separate  the  two  paths  vertically  as  in  the  figure  and 
assume  the  light  travelling  through  the  slab  along  the  upper  path  takes  longer 
than  for  a  ray  passing  through  along  the  lower  path.  We  denote  the  transit  time 
difference  by  A  Sr.  Due  to  the  transit  time  difference,  light  on  the  lower  path 
will  be  able  to  travel  slightly  beyond  the  end  of  the  slab  in  the  same  time  it  takes 
light  along  the  upper  path  to  transit  the  slab.  The  mean  additional  propagation 
distance  for  this  portion  of  the  beam  will  be  approximately  cASt/uq,  where  the 
contribution  due  to  V\  is  assumed  negligible  (n\  <C  no). 


Figure  4.  Ray  deflection  due  to  differential  time  to  cross  a  turbulent  slab. 


Due  to  the  additional  propagation  of  the  light  along  the  lower  path,  the 
orientation  of  the  wavefront  of  the  beam  will  be  slightly  deflected  by  amount 
Aft  upon  transit  through  the  slab.  This  assumes  that  the  separation  distance 
between  the  two  paths  (Ay)  is  small  enough  that  the  tilt  effect  is  the  primary 
modifying  effect  of  the  turbulence  between  these  two  points.  The  deflection 
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angle  can  then  be  computed  as  the  ratio  of  the  additional  propagation  length 
cASt/uq  to  the  separation  distance  dy: 


A  a  = 


e  A St 
n0  Ay  ‘ 


(129) 


Let  us  now  take  the  limit  of  this  procedure  as  Ay  approaches  zero.  Deflection 
A  a  remains  finite  upon  passage  of  the  ray  through  the  slab,  but  Ay  — >  dy  and 
A  St  —>  d.Sr  become  differentials.  The  differential  form  of  this  result  can  then 
be  generalized  to  a  two-dimensional  deflection,  such  that 

A  a  =  —  Vjl  St,  (130) 

n0 

where  the  gradient  operates  in  the  transverse  plane  perpendicular  to  the  main 
axis  of  propagation. 

Representing  St  in  terms  of  its  Fourier  series  coefficients  Tim ,  we  see  that  St 
can  be  evaluated  at  an  arbitrary  location  via 


OO  OC 

Sr(x,y)  =  E  E  Tim  exp 

i  =  — OO  771  — —  OC 


(131) 


In  this  form,  the  gradient  operator  can  be  evaluated  directly  within  the  Fourier 
representation  itself  to  produce 


A  a(x,y)  = 
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Of  course,  computationally  speaking,  we  can  only  evaluate  Aa  at  a  finite  series 
of  points.  Under  these  circumstances  an  FFT  can  be  used,  and  intermediate 
points  can  be  evaluated  using  interpolation  techniques.  Equation  (132)  provides 
a  means  for  evaluating  this  finite  grid  of  deflection  vectors. 

This  result  apparently  concludes  our  analysis.  However,  there  is  a  problem. 
The  gradient  operation  only  applies  to  infinitessimally  narrow  rays  of  light.  But 
the  light  received  at  a  single  pixel  site  within  the  image  plane  must  arise  from  a 
finite  area  within  the  object  plane.  Also,  the  energy  from  this  source  region  must 
pass  through  a  finite-sized  receiver  aperture.  These  factors  indicate  that  finite 
width  beams  must  be  considered  in  evaluating  Aa.  Elements  of  this  issue  were 
considered  previously  when  evaluating  beam-wave  angle-of-arrival  calculations 
in  Tofsted  (1992).  That  document  contained  a  derivation  for  angle-of-arrival 
variance  of  a  uniform  cross-section  beam.  In  that  analysis  the  deflection  angle 
was  determined  by  considering  the  tilt  in  the  beam  across  a  distance  D ,  which 
was  nominally  assigned  to  the  value  of  the  beam  diameter. 


A  similar  approach  can  be  used  here,  but  must  be  modified  in  light  of  the 
imaging  scenario  at  hand.  There  are  three  aspects  to  these  modifications: 
First,  in  equation  (132)  all  frequencies  of  distortions  were  included  in  the 
spectral  summation.  However,  for  turbulence  fluctuations  whose  wavelengths 
are  comparable  to  or  narrower  than  the  beam  width,  their  primary  effect  will 
be  to  distort  the  beam  structure  rather  than  to  tilt  its  direction  of  propagation. 
Even  wavelengths  larger  than  the  beam  width  will  be  less  effective  at  causing 
tilt  because  the  beam  width  is  a  sizable  fraction  of  the  wavelength.  Second,  we 
must  somehow  define  what  we  mean  when  referring  to  the  beam  width.  The 
emissions  from  the  object  surfaces  radiate  as  spherical  waves  from  a  multitude 
of  independent  (incoherent)  regions  within  the  footprint  of  a  single  pixel’s  IFOV 
(Instantaneous  Field  Of  View),  (f>p .  In  the  absence  of  excessive  turbulence 
scattering  of  the  propagated  wave,  we  may  designate  a  volume  that  defines 
the  space  through  which  the  energy  arising  from  the  pixel  region  flows  through 
to  enter  the  receiver  aperture.  It  is  the  characteristic  width  of  this  volume  that 
is  needed  to  describe  the  influence  of  a  particular  spatial  frequency  on  the  tilt 
factor  Aa.  Lastly,  we  may  inquire  regarding  the  shape  of  the  “beam”  as  it  travels 
through  space.  Since  it  is  not  uniformly  square  or  circular  (it  metamorphoses 
from  one  shape  to  the  other  over  the  course  of  propagating  from  the  object  plane 
to  the  receiver  aperture),  it  must  have  different  effects  for  every  different  spatial 
frequency  at  every  different  position  along  the  optical  axis.  Moreover,  even  if  we 
considered  the  initial  radiance  pattern  arising  from  a  particular  source  region 
in  the  object  plane  as  constant  in  magnitude,  because  of  the  morphological 
changes,  the  irradiance  pattern  would  not  remain  constant  within  its  envelope 
as  a  function  of  path  position. 

Under  ideal  conditions  we  would  assume  the  receptor  area  for  a  single  pixel 
would  be  some  square  region  on  the  receiver  plane  (assuming  square  pixels). 
Ignoring  system  blur  and  turbulence  effects,  we  would  project  that  the  energy 
falling  within  a  single  pixel  would  arise  from  an  area  <f>pZ2  on  the  object  plane. 
Assuming  each  point  on  the  object  plane  associated  with  this  pixel  emits  a 
spherical  wave  of  radiance  7o,  and  that  only  the  portion  of  this  wave  that  enters 
the  system  aperture  is  significant,  we  can  model  the  irradiance  pattern  of  a  pixel 
region  in  the  object  plane  as 
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M  <  D/ 2; 
|r|  >  D/2. 


(134) 


For  convenience  the  pixel  position  was  assumed  to  lie  along  the  optical  axis. 


The  meaning  of  this  function  is  that  a  scaled  version  of  the  original  projected 
pixel  shape  is  convolved  with  a  scaled  version  of  the  entrance  aperture  shape 
of  the  receiver.  We  can  parameterize  these  results  if  we  set  P  =  (ppZ  and 
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let  PQ  =  (1  -  z/Z )  P  be  the  scaled  width  of  a  square  at  distance  2  just  as 
Dq  =  (z/Z)  D  is  a  scaled  diameter  of  a  circle  at  the  same  position.  Using  this 
method,  we  have  an  irradiance  pattern  which  transitions  smoothly  between  a 
square  pattern  at  the  source  and  a  circular  pattern  at  the  receiver. 

We  can  now  create  two  dimensionless  variables  A  =  (Dq  ~~  Pq)/{Dq  A-  Pq)i 
B  =  (Dq  +  Pq)/Xt,  where  Ar  is  a  single  wavelength  of  turbulent  fluctuations 
(either  X/Z  or  Y/rri  in  equation  (132)).  Due  to  the  separability  of  the  turbulent 
effects  in  x  and  y.  we  can  evaluate  tilt  efficiencies  for  components  on  each  axis 
individually.  Similar  to  the  approach  used  by  Fried  (1965)  we  can  evaluate  an 
overall  tilt  effect  on  the  aggregate  of  irradiance  passing  through  a  given  layer. 
We  evaluate  this  quantity  as  follows:  From  equation  (132)  we  extract  the  (l,  m) 
and  (— Z,  —  m)  components  of  the  summation  and  combine  them  to  produce  a  tilt 
effect  for  a  single  characteristic  frequency.  Let  us  then  write  T)m  as  a  magnitude 
M  and  a  complex  phase  9 ,  T)m  =  M  etG .  Considering  only  the  x  axis  effects  we 
have  the  net  component  of  A5  in  the  x  direction  of 
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We  simplify  this  result  by  considering  only  a  component  along  the  x  axis  (m  — 
0).  This  step  is  necessary  to  simplify  the  mathematics,  but  can  be  seen  to  be 
possible  for  any  value  of  m  by  a  rotation  of  axes.  We  can  then  write 
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A  dimensional  analysis  of  equation  (132)  reveals  M  has  units  of  length.  To 
determine  the  mean  angular  tilt  of  a  beam,  wc  need  to  weigh  the  amount  of  the 
tilt  (here  Aaxim)  by  the  amount  of  energy  being  affected  at  every  point  in  the 
beam: 


A — (  \  lm(s)  Wh  z)  ds 

Aa(XT,  9)-  r  r  J(5  z)  dg 


(137) 


This  evaluation  determines  the  average  deflection  of  the  beam  for  a  particular 
value  of  0,  but  to  determine  the  overall  efficiency  of  the  turbulence  wavelength 
A t  in  deflecting  radiation,  this  result  needs  to  be  averaged  over  all  possible 
values  of  9.  Further,  since  A  a  is  a  zero  mean  variable,  it  will  be  necessary  to 
consider  the  second  moment  in  seeking  a  meaningful  statistic: 


2n 

^Aa2(Ay)^  =  —  ^  A  a  (Ay,  9)d9.  (138) 
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The  last  step  in  this  analysis  is  to  normalize  equation  (138)  by  dividing  by  the 
value  of  the  said  measure  when  the  irradiance  profile  approaches  zero  width. 
Evaluation  of  this  case  yields  the  limiting  value  M2/(2Xj<). 


However,  we  are  not  interested  in  the  efficiency  of  generating  A  a2,  but  A  a,  so 
the  final  efficiency  metric  chosen  should  be  proportional  to  the  square  root  of 


the  former  quantity: 


(139) 


Here,  efficiency  c  depends  on  both  A  and  B,  where  A  characterizes  the  shape 
of  the  irradiance  pattern  and  B  indicates  the  relative  size  of  the  pattern  to  the 
turbulence  wavelength.  For  the  case  of  a  square  pattern  A  =  —  1,  and  one  can 
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For  a  circular  pattern  j4  =  +1,  and  one  can  derive 


e  = 


7 xB 


(141) 


These  functions  are  very  similar  in  appearance  though  the  zeros  occur  at  different 
intervals  and  the  patterns  extinguish  at  different  rates. 

In  evaluating  the  more  general  cases,  wc  note  that  the  factor  M/Xt  drops  out 
when  dividing  by  1/2  M2/ Ay.  Second,  using  the  rule  sin(a  +  6)  =  sin(a)  cos(ft)  + 
cos  (a)  sin(6),  we  may  expand  the  term  sin  ((9  +  2nx/Xj-)  and  represent  the 
evaluation  in  terms  of  two  integrals  to  evaluate  equation  (137).  Of  these  two, 
one  will  contain  the  term  siii(27rx/Av),  and  because  sine  is  an  odd  function, 
while  the  irradiance  pattern  is  even,  the  net  integral  must  evaluate  to  zero. 

The  remaining  term  in  equation  (137)  involves  a  sin(0)  quantity  which  factors 
out  of  the  integral.  Equation  (138)  is  thus  considerably  simplified  and  yields  the 
result  (27t)_1  sin2(0)  dQ  =  1/2.  (Note:  This  1/2  factor  cancels  with  the  1/2 

in  the  normalizing  factor.)  We  then  further  normalize  the  integration  procedure 
by  the  changes  of  variables  x  =  XtBu  and  y  =  X  j-Bv.  These  substitutions  result 
in  an  argument  of  2ttBu  in  the  cosine  variable  and  to  normalized  limits  for  the 
integration  region  of  the  irradiance  signature: 

1/2 

e(A,  B)  =  y*  cos(2n Bu)  In(A,  B,  u)du  ,  (142) 

-i/2 


where  In  is  a  normalized  irradiance  pattern  involving  a  final  simplification 
whereby  the  v  variable  is  integrated  out.  The  resulting  integral  equation  can  be 
rapidly  evaluated  numerically  for  a  series  of  patterns  with  different  A  parameter 
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Normalized  Pattern  Position  (11) 


Figure  5.  Normalized  irradiance  profiles  plotted  for  various  values  of  the  parameter 
A.  The  area  under  each  curve  equals  unity.  The  normalized  profile  position  variable 
u  ranges  over  the  width  of  the  pattern  from  -1/2  to  +1/2. 


values.  Wc  present  a  scries  of  In  patterns  plotted  for  A  =  in  steps  of  1/4 

in  figure  5. 

Use  of  the  normalized  irradiance  profiles  permits  us  to  evaluate  the  efficiency 
factor  using  the  integral  relation 


e(A,  B)  = 


1/2 
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cos(27 xBu)  In(A,  u )  du 


(143) 


Resulting  efficiency  curves  for  the  cases  included  in  figure  5  are  plotted  in 
figure  6. 

Since  we  have  previously  noted  that  the  effects  of  beam  shape  can  be  assessed 
independently  for  each  axis,  we  now  write  the  final  equation  for  Ac*  by  inserting 
efficiency  factors  e(A,  B)  for  each  axis  into  our  A <3  equation: 


A  a(x,y)  =  —  EE  i  2n 
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where  C(z)  =■  PQ(z)  +  DQ(z)  is  introduced. 

Though  the  deflector  screen  model  cannot  handle  scintillation  effects,  influences 
of  wavefront  degradation  can  be  approximated  by  subsampling  the  arriving 
radiance  over  the  system  entrance  pupil.  This  subsampling  procedure  is 
necessary  when  the  turbulent  coherence  diameter  is  markedly  smaller  than  the 
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Figure  6.  Efficiency  response  function,  e.  plotted  as  a  function  of  the  shape  parameter 
A  and  the  pattern  ratio  B.  The  range  of  A  values  plotted  indicates  the  behavior  of 
the  efficiency  across  intermediate  conditions  spanning  between  the  limiting  cases. 

system  entrance  pupil  diameter.  In  effect,  different  portions  of  the  entrance 
pupil  are  seeing  light  arriving  from  a  single  direction  that  originated  from 
different  pixel  regions  in  the  object  plane.  Thus,  it  is  possible  to  include 
simulated  blurring  effects  by  subsampling  the  entrance  pupil.  To  determine  if 
such  subsampling  is  necessary,  a  coherence  diameter  (r0)  can  be  computed  and 
compared  with  the  system  entrance  pupil  diameter  ( D ).  If  the  ratio  D/r0  >  1, 
subsampling  will  be  necessary.  A  similar  comparison  can  be  made  between  a 
coherence  diameter  (r00),  calculated  using  a  path  weighting  function  oriented 
toward  the  object  plane,  and  the  width  P  =  Z(j>v .  If  P/r00  >  1,  it  will  be 
necessary  to  subsample  angular  trajectories  within  the  field  of  regard  of  each 
pixel. 

The  effect  of  subsampling  either  or  both  of  these  regions  would  be  to  reduce  P 
and/or  D  by  factors  y/N ~  and/or  y/Nd,  respectively.  Call  these  new  subsampled 
scales  Pf  and  D When  subsampling  the  receiver  aperture  it  would  be  necessary 
to  average  the  sample  results  obtained  to  produce  a  mean  result  for  each  pixel 
based  on  the  source  regions  mapped  to  by  the  iVd  samples.  From  the  above 
discussion,  it  would  be  necessary  to  consider  the  shape  of  the  subsamples  as 
well  as  their  number.  For  subsamples  in  the  source  pixel  region,  it  would  be 
more  advantageous  to  subdivide  it  into  equal-shaped  squares.  For  subsamples 
in  a  circular  receiver  aperture  it  might  be  more  advantageous  to  model  the  effects 
in  terms  of  hexagons  which  have  the  approximate  shape  of  circles. 

Having  performed  this  analysis,  we  may  be  able  to  determine  which  spectral 
region  has  the  most  effect  on  tilt.  For  the  moment  let  us  consider  the  case  of  small 
pixels  and  a  large  system  aperture.  In  this  case  the  shape  of  the  spectral  profile 
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fits  the  case  A  =  1  over  nearly  the  entire  path.  Hence,  based  on  equation  (144), 
the  magnitude  of  the  Aa  kernel  is  approximately  proportional  to 


|N|  OC  a  yJ^nW)  e[l,  crC(z)].  (145) 

In  this  equation  we  have  approximated  the  effect  of  either  of  the  two  orthogonal 
Bessel  functions  by  a  single  radially  symmetric  e  function.  We  have  also 
assessed  the  radial  behavior  of  the  gradient  operator  as  being  proportional  to 
the  magnitude  of  the  vector.  The  resulting  magnitude  is  just  the  modulus  of 
the  Fourier  transform  of  Aa,  here  written  as  R.  We  show  sample  plots  of  this 
function  in  figure  7.  We  have  plotted  curves  for  values  of  L0  that  arc  10 
100 Dz,  and  1000 Dz. 


Figure  7.  Proportional  magnitude  of  the  Fourier  spectrum  of  the  beam  deflections 
normalized  with  respect  to  the  outer  scale  length. 


As  can  be  seen,  finite  pixels  and  entrance  apertures  tend  to  negate  the  efficiency 
of  higher  frequency  turbulence  in  creating  beam  tilt.  However,  the  highest 
frequency  components  obviously  do  not  generate  the  greatest  influence  on  the 
tilt.  Rather,  it  is  spectral  components  in  the  vicinity  of  the  outer  scale  that 
dominate  tilt  effects. 

4.5  Method  Comparison 
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Before  concluding,  it  is  worthwhile  to  mention  that  the  temporal  cross 
correlation,  r*T,  derived  in  section  4.2  is  equivalent  to  the  optical  path  length 
cross  correlation  derived  by  Fleck  et  al.  (1975)  and  Martin  and  Flatte  (1988) 


except  for  changes  in  nomenclature.  Nevertheless,  there  is  a  slight  discontinuity 
between  these  results  and  those  of  Goodman  (1985). 

To  study  this  discontinuity  we  must  begin  by  translating  between  our 
nomenclature  and  that  used  by  these  other  authors.  Martin  and  Flatte  (1988) 
develop  a  phase  correlation  function,  Bg,  which  is  very  similar  to  the  F  (optical 
depth)  correlation  analysis  of  Fleck  et  al.  (1976).  The  use  of  the  B  convention 
for  covariances  is  standard  throughout  a  significant  portion  of  the  available 
literature  (cf.  Beland,  1993).  Martin  and  Flatte’s  Bg  is  proportional  to  our 
rsr  function: 

Bg(s)  =  k 2  c2  rgT(s).  (146) 

Bg  thus  equals  a  set  of  scaling  constants  times  the  2-D  inverse  transform  of 
$n(«t>  0).  Hence,  taking  the  forward  transform  of  Bg  produces  their  function 

^e(«x)  =  27rA;2  A$„(kj.,  0),  (147) 

where  iz±  is  the  transverse  plane  component  of  their  k  frequency  variable.  is 
thus  the  2-D  power  spectrum  for  phase  fluctuations. 

Previously  we  discussed  the  factor  2ttA ,  indicating  that  the  range  step  A  was 
bounded  both  above  and  below.  We  now  note  a  comparison  between  this  result 
and  a  result  for  a  similar  analysis  by  Goodman  (1985).  Goodman’s  power 
spectrum  of  phase  (Goodman’s  equation  (8.6-25))  is  written  as, 
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The  relationship  between  and  Fs  should  be  one-to-one,  but  Goodman’s  Fjj 
function  contains  a  term  that  depends  on  the  dimensionless  variable  it  =  K2±A/k. 
In  the  limit  as  u  — »  0  this  term  approaches  2,  which  agrees  with  Martin  and 
Flatte.  But  for  large  u.  the  term  fluctuates  about  unity.  Hence,  for  accurate 
evaluation ^)f  the  phase  screen  and  deflector  screen  approaches  we  must  ensure 
that  A  <C  k/K2±  in  addition  to  our  other  conditions.  Goodman’s  analysis  is  also 
instructive  in  that  his  results  indicate  that  the  mean  wavenumber  of  the  receiver 
has  an  influence  on  the  overall  requirements  of  the  wave  propagation  methods 
used. 
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5.  Conclusions 


In  this  document  two  key  methodologies  for  simulating  the  propagation  of 
imagery  through  optical  turbulence  have  been  presented  in  a  comprehensive 
analysis.  The  mathematical  basis  for  such  analysis  rests  upon  the  various 
forms  of  the  Fourier  transform,  including  Fourier  series  expansions  for  periodic 
functions  and  FFTs  for  sampled  functions.  Based  on  these  Fourier  methods 
it  was  shown  that  various  means  of  representing  the  turbulent  refractive  index 
power  spectrum  exist,  and  the  interrelationships  between  these  expressions  were 
described,  particularly  between  the  form  commonly  used  in  turbulence  analysis 
literature  and  the  form  useful  in  FFT  calculations. 

A  derivation  was  also  presented  for  evaluating  the  propagation  statistics  of 
cross  correlation  of  fluctuations  in  the  time  of  transit  of  wave  energy  through 
a  turbulent  layer  of  thickness  A.  Two  methods  were  described,  one  supporting 
wave  optical  calculations  known  as  the  phase  screen  approach,  and  a  second 
method  based  on  raytracing  called  the  deflector  screen  approach. 

The  phase  screen  method  is  more  accurate  because  it  treats  the  propagating 
energy  using  wave  optics.  However,  this  method  may  be  prohibitively  more 
expensive  in  computational  time  as  the  fluctuations  in  the  field  must  be  tracked 
with  high  precision  at  each  integration  step.  By  comparison,  the  tilt  effects 
modeled  by  the  deflector  screen  method  account  for  the  largest  energy  portion 
of  the  turbulence  spectrum  and  should  be  sufficient  to  capture  large  scale 
distortions  of  objects  due  to  turbulence.  The  limitation  is  that  the  deflector 
screen  method  cannot  properly  handle  turbulent  scintillation  or  image  blurring 
effects  because  only  the  tilt  of  the  beam  is  considered,  not  turbulent  destruction 
of  the  propagating  wavefront,  or  loss  of  wavefront  coherence.  For  far-IR,  systems 
such  effects  arc  not  critical.  To  determine  the  limits  of  applicability  of  the 
deflector  screen  method  it  is  recommended  that  the  phase  screen  approach  be 
used  as  a  benchmark.  With  regard  to  processing  requirements  the  deflector 
screen  method  will  likely  provide  sufficient  accuracy  while  providing  rapid 
processing  capability  in  support  of  scenario  generation  for  user  perception 
testing.  The  phase  screen  method,  on  the  other  hand,  will  require  very  high- 
powered  computer  capabilities  in  support  of  realistic  propagation  studies.  The 
result  will  likely  be  that  the  phase  screen  method  will  see  limited  use,  but  will 
provide  valuable  information,  while  the  bulk  of  the  processing  is  accomplished 
using  the  deflector  screen  method  which  is  computationally  much  more  efficient. 
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Acronyms 


EOSAEL 

Electro-Optical  Systems  Atmospheric  Effects  Library 

FLIR 

Forward-Looking  Infrared 

FFT 

Fast  Fourier  Transform 

IFOV 

instantaneous  field  of  view 

IR 

infrared 

LSI 

linear  shift  invariant 

MFT 

modulation  transfer  function 

RHS 

right  hand  side 

55 


Distribution 


Copies 

NASA  MARSHALL  SPACE  FLT  CTR  1 

ATMOSPHERIC  SCIENCES  DIV  SDO 1 
ATTN  DR  FICHTL 
HUNTSVILLE  AL  35802 

NASA  SPACE  FLT  CTR  1 

ATMOSPHERIC  SCIENCES  DIV 
CODE  ED  41  1 
HUNTSVILLE  AL  35812 

US  ARMY  MISSILE  CMND  1 

AMSMI  RD  AS  SS 

ATTN  MR  H  F  ANDERSON 

REDSTONE  ARSENAL  AL  35898-5253 

US  ARMY  MISSILE  CMND  1 

AMSMI  RD  AS  SS 

ATTN  MR  B  WILLIAMS 

REDSTONE  ARSENAL  AL  35898-5253 

US  ARMY  MISSILE  CMND  1 

AMSMI  RD  DE  SE 

ATTN  MR  GORDON  LILL  JR 

REDSTONE  ARSENAL  AL  35898-5245 

US  ARMY  MISSILE  CMND  1 

REDSTONE  SCI  INFO  CTR 

AMSMI  RD  CS  R  DOC 

REDSTONE  ARSENAL  AL  35898-5241 

US  ARMY  MISSILE  CMND  1 

AMSMI 

REDSTONE  ARSENAL  AL  35898-5253 


PACIFIC  MISSILE  TEST  CTR 
GEOPHYSICS  DIV 
ATTN  CODE  3250 
POINT  MUGU  CA  93042-5000 

ATMOSPHERIC  PROPAGATION  BRANCH 
SPAWARSYSCEN  SAN  DIEGO  D858 
49170  PROPAGATION  PATH 
SAN  DIEGO  CA  92152-7385 

METEOROLOGIST  IN  CHARGE 
KWAJALEIN  MISSILE  RANGE 
PO  BOX  67 

APO  SAN  FRANCISCO  CA  96555 

NCAR  LIBRARY  SERIALS 
NATL  CTR  FOR  ATMOS  RSCH 
PO  BOX  3000 
BOULDER  CO  80307-3000 

HEADQUARTERS  DEPT  OF  ARMY 
DAMI  POI 
ATTN  LEE  PAGE 
WASHINGTON  DC  20310-1067 

US  ARMY  INFANTRY 
ATSH  CD  CS  OR 
ATTN  DR  E  DUTOIT 
FT  BENNING  GA  30905-5090 

HQ  AFWA/DNX 

106  PEACEKEEPER  DR  STE  2N3 
OFFUTTAFB  NE  681 13-4039 

PHILLIPS  LABORATORY 
PLLYP 

ATTN  MR  CHISHOLM 
HANSCOM  AFB  MA  01731-5000 

PHILLIPS  LABORATORY 
PL  LYP  3 

HANSCOM  AFB  MA  01731-5000 

AFRL/VSBL 
29  RANDOLPH  RD 
HANSCOM  AFB  MA  01731 


58 


ARL  CHEMICAL  BIOLOGY 
NUC  EFFECTS  DIV 
AMSRL  SL  CO 
APG  MD  21010-5423 

US  ARMY  MATERIEL  SYST 
ANALYSIS  ACTIVITY 
AMSXY 

APG  MD  21005-5071 

US  ARMY  RESEARCH  LABORATORY 
AMSRL  D 

2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1145 

US  ARMY  RESEARCH  LABORATORY 
AMSRL  OP  Cl  SD  TL 
2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1145 

US  ARMY  RESEARCH  LABORATORY 

AMSRL  Cl  LL 

ADELPHI  MD  20783-1197 

US  ARMY  RESEARCH  LABORATORY 
AMSRL  SS  SH 
ATTN  DR  SZTANKAY 
2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1 145 

US  ARMY  RESEARCH  LABORATORY 

AMSRL  Cl 

ATTN  J  GANTT 

2800  POWDER  MILL  ROAD 

ADELPHI  MD  20783-1 197 

US  ARMY  RESEARCH  LABORATORY 
AMSRL 

2800  POWDER  MILL  ROAD 
ADELPHI  MD  20783-1145 

NATIONAL  SECURITY  AGCY  W21 
ATTN  DR  LONGBOTHUM 
9800  SAVAGE  ROAD 
FT  GEORGE  G  MEADE  MD  20755-6000 


US  ARMY  RSRC  OFC 
ATTN  AMXRO  GS  DR  BACH 
PO  BOX  12211 
RTP  NC  27009 

DR  JERRY  DAVIS 
NCSU 

PO  BOX  8208 
RALEIGH  NC  27650-8208 

US  ARMY  CECRL 
CECRL  GP 
ATTN  DR  DETSCH 
HANOVER  NH  03755-1290 

US  ARMY  ARDEC 
SMCAR  IMI I  BLDG  59 
DOVER  NJ  07806-5000 

ARMY  DUGWAY  PROVING  GRD 
STEDP  MT  DA  L  3 
DUGWAY  UT  84022-5000 

ARMY  DUGWAY  PROVING  GRD 
STEDP  MT  M 
ATTN  MR  BOWERS 
DUGWAY  UT  84022-5000 

DEPT  OF  THE  AIR  FORCE 

OL  A  2D  WEATHER  SQUAD  MAC 

HOLLOMAN  AFB  NM  88330-5000 

PL  WE 

KIRTLAND  AFB  NM  87118-6008 

USAF  ROME  LAB  TECH 
CORRIDOR  W  STE  262  RL  SUL 
26  ELECTR  PKWY  BLD  106 
GRIFFISS  AFB  NY  13441-4514 

AFMC  DOW 

WRIGHT  PATTERSON  AFB  OH  45433-5000 


60 


US  ARMY  FIELD  ARTILLERY  SCHOOL  1 

ATSFTSMTA 

FT  SILL  OK  73503-5600 

US  ARMY  FOREIGN  SCI  TECH  CTR  1 

CM 

220  7TH  STREET  NE 
CHARLOTTESVILLE  VA  22902-5396 

NAVAL  SURFACE  WEAPONS  CTR  1 

CODE  G63 

DAHLGREN  VA  22448-5000 

US  ARMY  OEC  1 

CSTE  EFS 

PARK  CENTER  IV 

4501  FORD  AVE 

ALEXANDRIA  VA  22302-1458 

US  ARMY  CORPS  OF  ENGRS  1 

ENGR  TOPOGRAPHICS  LAB 
ETL  GS  LB 

FT  BELVOIR  VA  22060 


US  ARMY  TOPO  ENGR  CTR  1 

CETEC  ZC  1 

FT  BELVOIR  VA  22060-5546 

SCI  AND  TECHNOLOGY  1 

101  RESEARCH  DRIVE 
HAMPTON  VA  23666-1340 

US  ARMY  NUCLEAR  CML  AGCY  1 

MONA  ZB  BLDG  2073 
SPRINGFIELD  VA  22150-3198 

USATRADOC  1 

ATCD  FA 

FT  MONROE  VA  23651-5170 

ATRC  WSS  R  1 

WSMR  NM  88002-5502 


) 


61 


US  ARMY  RESEARCH  LABORATORY 
AMSRL  Cl  E 
COMP  &  INFO  SCI  DIR 
WSMR  NM  88002-5501 

DTIC 

8725  JOHN  J  KINGMAN  RD 
STE  0944 

FT  BELVOIR  VA  22060-6218 

US  ARMY  MISSILE  CMND 
AMSMI 

REDSTONE  ARSENAL  AL  35898-5243 

US  ARMY  DUGWAY  PROVING  GRD 
STEDP3 

DUGWAY  UT  84022-5000 

USTRADOC 
ATCD  FA 

FT  MONROE  VA  23651-5170 

WSMR  TECH  LIBRARY  BR 
STEWS  IM  IT 
WSMR  NM  88002 

US  ARMY  RESEARCH  LAB 
AMSRL  D  DR  D  SMITH 
2800  POWDER  MILL  RD 
ADELPHI  MD  20783-1 197 

US  ARMY  CECOM 
INFORMATION  &  INTELLIGENCE 
WARFARE  DIRECTORATE 
ATTN  AMSEL  RD  IW  IP 
FORT  MONMOUTH  NJ  07703-521 1 

US  ARMY  RESEARCH  LAB 
ATTN  AMSRL  Cl  EW 
MR  TOFSTED 
WSMR  NM  88002-5513 


Record  copy 

TOTAL 


62 


